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Abstract 

Background:  The  hypothalamic-pituitary-adrenal  (HPA)  axis  is  a  neuroendocrine  system  that  regulates  numerous 
physiological  processes.  Disruptions  in  the  activity  of  the  HPA  axis  are  correlated  with  stress-related  diseases  such  as 
post-traumatic  stress  disorder  (PTSD)  and  major  depressive  disorder.  In  this  paper,  we  characterize  "normal"  and 
"diseased"  states  of  the  HPA  axis  as  basins  of  attraction  of  a  dynamical  system  describing  the  inhibition  of  peptide 
hormones  such  as  corticotropin-releasing  hormone  (CRH)  and  adrenocorticotropic  hormone  (ACTH)  by  circulating 
glucocorticoids  such  as  cortisol  (CORT). 

Results:  In  addition  to  including  key  physiological  features  such  as  ultradian  oscillations  in  cortisol  levels  and 
self-upregulation  of  CRH  neuron  activity,  our  model  distinguishes  the  relatively  slow  process  of  cortisol-mediated  CRH 
biosynthesis  from  rapid  trans-synaptic  effects  that  regulate  the  CRH  secretion  process.  We  show  that  the  slow 
component  of  the  negative  feedback  allows  external  stress-induced  reversible  transitions  between  "normal"  and 
"diseased"  states  in  novel  intensity-,  duration-,  and  timing-dependent  ways. 

Conclusion:  Our  two-step  negative  feedback  model  suggests  a  mechanism  whereby  exposure  therapy  of  stress 
disorders  such  as  PTSD  may  act  to  normalize  downstream  dysregulation  of  the  HPA  axis.  Our  analysis  provides  a 
causative  rationale  for  improving  treatments  and  guiding  the  design  of  new  protocols. 
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Background 

Stress  is  an  essential  component  of  an  organisms  attempt 
to  adjust  its  internal  state  in  response  to  environmen¬ 
tal  change.  The  experience,  or  even  the  perception  of 
physical  and/or  environmental  change,  induces  stress 
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responses  such  as  the  secretion  of  glucocorticoids  hor¬ 
mones  (CORT)  -  cortisol  in  humans  and  corticosterone 
in  rodents  -  by  the  adrenal  gland.  The  adrenal  gland 
is  one  component  of  the  hypothalamic-pituitary-adrenal 
(HPA)  axis,  a  collection  of  interacting  neuroendocrine 
cells  and  endocrine  glands  that  play  a  central  role  in 
stress  response.  The  basic  interactions  involving  the  HPA 
axis  are  shown  in  Fig.  1.  The  paraventricular  nucleus 
(PVN)  of  the  hypothalamus  receives  synaptic  inputs  from 
various  neural  pathways  via  the  central  nervous  system 
that  are  activated  by  both  cognitive  and  physical  stres¬ 
sors.  Once  stimulated,  CRH  neurons  in  the  PVN  secrete 
corticotropin-releasing  hormone  (CRH),  which  then  stim¬ 
ulates  the  anterior  pituitary  gland  to  release  adrenocorti- 
cotropin  hormone  (ACTH)  into  the  bloodstream.  ACTH 
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Fig.  1  Schematic  of  HPA  axis,  a  Stress  is  processed  in  the  central  nervous  system  (CNS)  and  a  signal  is  relayed  to  the  PVN  in  the  hypothalamus  to 
activate  CRH  secretion  into  the  hypophyseal  portal  system,  b  CRH  diffuses  to  the  pituitary  gland  and  activates  ACTH  secretion.  ACTH  travels  down  to 
the  adrenal  cortex  to  activate  cortisol  (CORT)  release.  Cortisol  inhibits  both  CRH  and  ACTH  secretion  to  down-regulate  its  own  production,  forming  a 
closed  loop.  In  the  pituitary  gland,  cortisol  binds  to  glucocorticoid  receptors  (GR)  (yellow box)  to  inhibit  ACTH  and  self-upregulate  GR  production. 
This  part  of  the  axis  comprises  the  PA  subsystem,  c  Negative  feedback  of  cortisol  affects  the  synthesis  process  in  the  hypothalamus,  which  indirectly 
suppresses  the  release  of  CRH.  External  inputs  such  as  stressors  and  circadian  inputs  also  directly  affect  the  release  rate  of  CRH 


then  activates  a  complex  signaling  cascade  in  the  adrenal 
cortex,  which  ultimately  releases  glucocorticoids  (Fig.  lb). 
In  return,  glucocorticoids  exert  a  negative  feedback  on 
the  hypothalamus  and  pituitary,  suppressing  CRH  and 
ACTH  release  and  synthesis  in  an  effort  to  return  them 
to  baseline  levels.  Classic  stress  responses  include  tran¬ 
sient  increases  in  levels  of  CRH,  ACTH,  and  cortisol. 
The  basic  components  and  organization  of  the  vertebrate 
neuroendocrine  stress  axis  arose  early  in  evolution  and 
the  HPA  axis,  in  particular,  has  been  conserved  across 
mammals  [1]. 

Dysregulation  in  the  HPA  axis  is  known  to  correlate 
with  a  number  of  stress-related  disorders.  Increased  corti¬ 
sol  (hypercortisolism)  is  associated  with  major  depressive 
disorder  (MDD)  [2,  3],  while  decreased  cortisol  (hypocor- 
tisolism)  is  a  feature  of  post-traumatic  stress  disorder 
(PTSD),  post  infectious  fatigue,  and  chronic  fatigue  syn¬ 
drome  (CFS)  [4-7].  Since  PTSD  develops  in  the  aftermath 
of  extreme  levels  of  stress  experienced  during  traumatic 
incidents  like  combat,  sexual  abuse,  or  life-threatening 
accidents,  its  progression  may  be  strongly  correlated  with 
disruption  of  the  HPA  axis  caused  by  stress  response.  For 
example,  lower  peak  and  nadir  cortisol  levels  were  found 
in  patients  with  combat- related  PTSD  [8]. 

Mathematical  models  of  the  HPA  axis  have  been  pre¬ 
viously  formulated  in  terms  of  dynamical  systems  of 
ordinary  differential  equations  (ODEs)  [9-12]  or  delay 
differential  equations  (DDEs)  [13-15]  that  describe  the 
time-evolution  of  the  key  regulating  hormones  of  the  HPA 
axis:  CRH,  ACTH,  and  cortisol.  These  models  [13,  14,  16] 
incorporate  positive  self-regulation  of  glucocorticoid 


receptor  expression  in  the  pituitary,  which  may  generate 
bistability  in  the  dynamical  structure  of  the  model  [17]. 
Of  the  two  stable  equilibrium  states,  one  is  characterized 
by  higher  levels  of  cortisol  and  is  identified  as  the  “nor¬ 
mal”  state.  The  other  is  characterized  by  lower  levels  of 
cortisol  and  can  be  interpreted  as  one  of  the  “diseased” 
states  associated  with  hypocortisolism.  Stresses  that  affect 
the  activity  of  neurons  in  the  PVN  are  described  as  pertur¬ 
bations  to  endogenous  CRH  secretion  activity.  Depending 
on  the  length  and  magnitude  of  the  stress  input,  the  sys¬ 
tem  may  or  may  not  shift  from  the  basin  of  attraction  of 
the  normal  steady  state  towards  that  of  the  diseased  one. 
If  such  a  transition  does  occur,  it  may  be  interpreted  as  the 
onset  of  disease.  A  later  model  [16]  describes  the  effect  of 
stress  on  the  HPA  axis  as  a  gradual  change  in  the  param¬ 
eter  values  representing  the  maximum  rate  of  CRH  pro¬ 
duction  and  the  strength  of  the  negative  feedback  activity 
of  cortisol.  In  this  model,  cortisol  secretion  patterns  are 
assumed  to  depend  solely  on  physiological  changes  arising 
from  e.g.,  anatomical  or  biochemical  changes  in  cells  or 
tissues.  Such  structural-level  variations  can  be  mathemat¬ 
ically  represented  by  changes  in  physiological  parameter 
values. 

These  two  classes  of  models  imply  qualitatively  different 
time  courses  of  disease  progression  [16,  17].  The  former 
suggests  that  the  abnormal  state  is  a  pre-existing  basin  of 
attraction  of  a  dynamical  model  that  stays  dormant  until  a 
sudden  transition  is  triggered  by  exposure  to  trauma  [17]. 
In  contrast,  the  latter  assumes  that  the  abnormal  state  is 
reached  by  the  slow  development  of  structural  changes  in 
physiology  due  to  the  traumatic  experience  [16].  Although 
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both  models  [16,  17]  describe  changes  in  hormonal  lev¬ 
els  experienced  by  PTSD  patients,  they  both  fail  to  exhibit 
stable  ultradian  oscillations  in  cortisol,  which  is  known  to 
play  a  role  in  determining  the  responsiveness  of  the  HPA 
axis  to  stressors  [18]. 

In  this  study,  we  consider  a  number  of  distinctive  physi¬ 
ological  features  of  the  HPA  axis  that  give  a  more  complete 
picture  of  the  dynamics  of  stress  disorders  and  that  have 
not  been  considered  in  previous  mathematical  models. 
These  include  the  effects  of  intrinsic  ultradian  oscilla¬ 
tions  on  HPA  dysregulation,  distinct  rapid  and  slow  feed¬ 
back  actions  of  cortisol,  and  the  correlation  between  HPA 
imbalance  and  disorders  induced  by  external  stress.  As 
with  the  majority  of  hormones  released  by  the  body,  corti¬ 
sol  levels  undergo  a  circadian  rhythm,  starting  low  during 
night  sleep,  rapidly  rising  and  reaching  its  peak  in  the  early 
morning,  then  gradually  falling  throughout  the  day.  Super¬ 
posed  on  this  slow  diurnal  cycle  is  an  ultradian  rhythm 
consisting  of  approximately  hourly  pulses.  CRH,  ACTH, 
and  cortisol  are  all  secreted  episodically,  with  the  pulses  of 
ACTH  slightly  preceding  those  of  cortisol  [19]. 

As  for  many  other  hormones  such  as  gonadotropin¬ 
releasing  hormone  (GnRH),  insulin,  and  growth  hormone 
(GH),  the  ultradian  release  pattern  of  glucocorticoids  is 
important  in  sustaining  normal  physiological  functions, 
such  as  regulating  gene  expression  in  the  hippocampus 
[20].  It  is  unclear  what  role  oscillations  play  in  home¬ 
ostasis,  but  the  time  of  onset  of  a  stressor  in  relation 
to  the  phase  of  the  ultradian  oscillation  has  been  shown 
to  influence  the  physiological  response  elicited  by  the 
stressor  [21]. 

To  distinguish  the  rapid  and  slow  actions  of  cortisol, 
we  separate  the  dynamics  of  biosynthesis  of  CRH  from 
its  secretion  process,  which  operate  over  very  different 
timescales  [22].  While  the  two  processes  are  mostly  inde¬ 
pendent  from  each  other,  the  rate  of  CRH  secretion  should 
depend  on  the  synthesis  process  since  CRH  peptides  must 
be  synthesized  first  before  being  released  (Fig.  lc).  On  the 
other  hand,  the  rate  of  CRH  peptide  synthesis  is  influ¬ 
enced  by  cortisol  levels,  which  in  turn,  are  regulated  by 
released  CRH  levels.  We  will  investigate  how  the  sep¬ 
aration  and  coupling  of  these  two  processes  can  allow 
stress-induced  dysregulations  of  the  HPA  axis. 

The  mathematical  model  we  derive  incorporates  the 
above  physiological  features  and  reflects  the  basic  physi¬ 
ology  of  the  HPA  axis  associated  with  delays  in  signaling, 
fast  and  slow  negative  feedback  mechanisms,  and  CRH 
self-upregulation  [23].  Within  an  appropriate  parameter 
regime,  our  model  exhibits  two  distinct  stable  oscillat¬ 
ing  states,  of  which  one  is  marked  by  a  larger  oscillation 
amplitude  and  a  higher  base  cortisol  level  than  the  other. 
These  two  states  will  be  referred  to  as  normal  and  dis¬ 
eased  states.  Our  interpretation  is  reminiscent  of  the 
two-state  dynamical  structure  that  arises  in  the  classic 


Fitzhugh-Nagumo  model  of  a  single  neuron,  in  which  rest¬ 
ing  and  spiking  states  emerge  as  bistable  modes  of  the 
model  [24],  or  in  models  of  neuronal  networks  where 
an  “epileptic  brain”  is  described  in  terms  of  the  distance 
between  a  normal  and  a  seizure  attractor  in  phase-space 
[25].  Our  main  result  is  that  stress-driven  transitions 
between  normal  and  diseased  states  can  arise  when  a  two- 
stage  negative  feedback  (of  cortisol  on  CRH)  mechanism 
is  incorporated.  The  possibility  of  such  transitions  lead  to 
a  number  of  novel  features  in  the  overall  system. 

Methods 

Models  of  HPA  dynamics  [13,  14,  16,  17,  26]  are  typi¬ 
cally  expressed  in  terms  of  ordinary  differential  equations 
(ODEs): 


%  =  pcI(T)fc(0)  -  dc(C), 
a  1 

(1) 

^ =  pACfA(OR,0)-dA(A ), 

(2) 

dO 

—  =poA(T )  -  do(0), 
al 

(3) 

d R 

-7^  =PrSr{OR)  -  di>(R), 
al 

(4) 

where  C(T),A(T),  and  O(T)  denote  the  plasma  concen¬ 
trations  of  CRH,  ACTH,  and  cortisol  at  time  T,  respec¬ 
tively.  R(T )  represents  the  availability  of  glucocorticoid 
receptor  (GR)  in  the  anterior  pituitary.  The  amount  of 
cortisol  bound  GR  is  typically  in  quasi-equilibrium  so 
concentration  of  the  ligand-receptor  complex  is  approx¬ 
imately  proportional  to  the  product  0(T)R(T)  [17].  The 
parameters  pa  (a  e  {C,A,  0,Rj)  relate  the  production 
rate  of  each  species  a  to  specific  factors  that  regulate 
the  rate  of  release/synthesis.  External  stresses  that  drive 
CRH  release  by  the  PVN  in  the  hypothalamus  are  rep¬ 
resented  by  the  input  signal  I(T).  The  function  /c(O) 
describes  the  negative  feedback  of  cortisol  on  CRH  levels 
in  the  PVN  while  /a  (OR,  O )  describes  the  negative  feed¬ 
back  of  cortisol  or  cortisol-GR  complex  (at  concentration 
0(T)R(T))  in  the  pituitary.  Both  are  mathematically  char¬ 
acterized  as  being  positive,  decreasing  functions  so  that 
Ia,c(‘)  >  0  and/^  c(’)  <  On  other  hand,  the  func¬ 
tion  gx(OR)  describes  the  self-upregulation  effect  of  the 
cortisol-GR  complex  on  GR  production  in  the  anterior 
pituitary  [27].  In  contrast  to  fA,c(‘)>  Sr(’)  is  a  positive  but 
increasing  function  of  OR  so  that  £/?(•)  >  0  andg^(-)  >  0. 
Finally,  the  degradation  functions  da(-)  describe  how  each 
hormone  and  receptor  is  cleared  and  may  be  linear  or 
nonlinear. 

Without  including  the  effects  of  the  glucocorticoid 
receptor  (neglecting  Eq.  4  and  assuming  /a( OR,  O)  = 
Ja(0)  in  Eq.  2),  Eqs.  1-3  form  a  rudimentary  “minimal” 
model  of  the  HPA  axis  [9,  28].  If/4, c(0  are  Hill-type  feed¬ 
back  functions  dependent  only  on  0(T)  and  da(-)  are 
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linear,  a  unique  global  stable  point  exists.  This  equilibrium 
point  transitions  to  a  limit  cycle  through  a  Hopf  bifurca¬ 
tion  but  only  within  nonphysiological  parameter  regimes 
[9].  The  inclusion  of  GR  and  its  self-upregulation  in 
the  anterior  pituitary  [17]  creates  two  stable  equilibrium 
states  of  the  system,  but  still  does  not  generate  oscillatory 
behavior.  More  recent  studies  extend  the  model  (rep¬ 
resented  by  Eqs.  1-4)  to  include  nonlinear  degradation 
[16]  or  constant  delay  to  account  for  delivery  of  ACTH 
and  synthesis  of  glucocorticoid  in  the  adrenal  gland  [13]. 
These  two  extended  models  exhibit  only  one  intrinsic  cir¬ 
cadian  [16]  or  ultradian  [13]  oscillating  cycle  for  any  given 
set  of  parameter  values,  precluding  the  interpretation  of 
normal  and  diseased  states  as  bistable  oscillating  modes  of 
the  model. 

Here,  we  develop  a  new  model  of  the  HPA  axis  by 
first  adapting  previous  work  [13]  where  a  physiologically- 
motivated  delay  was  introduced  into  Eq.  3,  giving  rise  to 
the  observed  ultradian  oscillations  [13].  We  then  improve 
the  model  by  distinguishing  the  relatively  slow  mecha¬ 
nism  underlying  the  cortisol-mediated  CRH  biosynthesis 
from  the  rapid  trans- synaptic  effects  that  regulate  CRH 
secretion.  This  allows  us  to  decompose  the  dynamics  into 
slow  and  fast  components.  Finally,  self-upregulation  of 
CRH  release  is  introduced  which  allows  for  bistability. 
These  ingredients  can  be  realistically  combined  in  a  way 
that  leads  to  novel,  clinically  identifiable  features  and  are 
systematically  developed  below. 

Ultradian  rhythm  and  time  delay 

Experiments  on  rats  show  a  3-6  min  inherent  delay  in 
the  response  of  the  adrenal  gland  to  ACTH  [29].  More¬ 
over,  in  experiments  performed  on  sheep  [30],  persistent 
ultradian  oscillations  were  observed  even  after  surgically 
removing  the  hypothalamus,  implying  that  oscillations  are 
inherent  to  the  pituitary-adrenal  (PA)  subsystem.  Since 
oscillations  can  be  induced  by  delays,  we  assume,  as 
in  Walker  et  al.  [13],  a  time  delay  T&  in  the  ACTH- 
mediated  activation  of  cortisol  production  downstream  of 
the  hypothalamus.  Equation  3  is  thus  modified  to 

^  =PoA(T-  Td)-d00.  (5) 

al 

Walker  et  al.  [13]  show  that  for  fixed  physiological  levels 
of  CRH,  the  solution  to  Eqs.  2,  4  and  5  leads  to  oscillatory 
A(T),  0(T),  and  R(T).  In  order  to  describe  the  observed 
periodic  cortisol  levels  in  normal  and  diseased  states,  the 
model  requires  two  oscillating  stable  states.  We  will  see 
that  dual  oscillating  states  can  arise  within  our  model 
when  the  delay  in  ACTH-mediated  activation  of  corti¬ 
sol  production  is  coupled  with  other  known  physiological 
processes  that  we  describe  below. 


Synthesis  of  CRH 

CRH  synthesis  involves  various  pathways,  including  CRH 
gene  transcription  and  transport  of  packaged  CRH  from 
the  cell  body  (soma)  to  their  axonal  terminals  where  they 
are  stored  prior  to  release.  Changes  in  the  steady  state 
of  the  synthesis  process  typically  occur  on  a  timescale  of 
minutes  to  hours.  On  the  other  hand,  the  secretory  release 
process  depends  on  changes  in  membrane  potential  at 
the  axonal  terminal  of  CRH  neurons,  which  occur  over 
millisecond  to  second  timescales. 

To  model  the  synthesis  and  release  process  separately, 
we  distinguish  two  compartments  of  CRH:  the  concentra¬ 
tion  of  stored  CRH  within  CRH  neurons  will  be  denoted 
CS(T),  while  levels  of  released  CRH  in  the  portal  vein  out¬ 
side  the  neurons  will  be  labeled  C(T)  (Fig.  lc).  Newly 
synthesized  CRH  will  first  be  stored,  thus  contributing 
to  Cs.  We  assume  that  the  stored  CRH  level  Cs  relaxes 
toward  a  target  value  set  by  the  function  Coo(O): 

dQ  =  0^(0)  -  Cs 
d  T  Tc  ' 

Here,  Tc  is  a  characteristic  time  constant  and  Coo(O) 
is  the  cortisol- dependent  target  level  of  stored  CRH. 
Equation  6  also  assumes  that  the  relatively  small  amounts 
of  CRH  released  into  the  bloodstream  do  not  significantly 
deplete  the  Cs  pool.  Note  that  the  effects  induced  by 
changing  cortisol  levels  are  immediate  as  the  production 
term  Coq(0)/Tc  is  adjusted  instantaneously  to  current 
cortisol  levels.  Our  model  thus  does  not  exclude  corti¬ 
sol  rapidly  acting  on  the  initial  transcription  activity,  as 
suggested  by  CRH  hnRNA  (precursor  mRNA)  measure¬ 
ments  [31].  On  the  other  hand,  the  time  required  to  reach 
the  steady  state  for  the  completely  synthesized  CRH  pep¬ 
tide  will  depend  on  the  characteristic  time  scale  constant 
Tc .  Ideally,  Tc  should  be  estimated  from  measurements 
of  the  pool  size  of  releasable  CRH  at  the  axonal  termi¬ 
nals.  To  best  of  our  knowledge,  there  are  currently  no 
such  measurements  available,  so  we  base  our  estimation 
on  mRNA  level  measurements.  We  believe  this  is  a  better 
representation  of  releasable  CRH  than  hnRNA  levels  since 
mRNA  synthesis  is  a  further  downstream  process.  Pre¬ 
vious  studies  have  shown  that  variations  in  CRH  mRNA 
due  to  changes  in  cortisol  levels  take  at  least  twelve  hours 
to  detect  [32].  Therefore,  we  estimate  Tc  >  12  hrs  = 
720  min.  The  negative  feedback  of  cortisol  on  CRH  lev¬ 
els  thus  acts  through  the  production  function  Coo(O)  on 
the  relatively  slow  timescale  Tc.  To  motivate  the  func¬ 
tional  form  of  Coo(O),  we  invoke  experiments  on  rats 
whose  adrenal  glands  had  been  surgically  removed  and 
in  which  glucocorticoid  levels  were  subsequently  kept 
fixed  (by  injecting  exogenous  glucocorticoid)  for  5-7  days 
[22,  33].  The  measured  CRH  mRNA  levels  in  the  PVN 
were  found  to  decrease  exponentially  with  the  level 
of  administered  glucocorticoid  [22,  33].  Assuming  the 


Ki  m  et  ai  Biology  Direct  (20 1 6)  1 1 : 1 3 


Page  5  of  18 


amount  of  releasable  CRH  is  proportional  to  the  amount 
of  measured  intracellular  CRH  mRNA,  we  can  approx¬ 
imate  Coo(O)  as  a  decreasing  exponential  function  of 
cortisol  level  O. 

Secretion  of  CRH 

To  describe  CRH  secretion,  we  consider  the  follow¬ 
ing  three  factors:  synaptic  inputs  to  CRH  cells  in  the 
PVN,  availability  of  releasable  CRH  peptide,  and  self- 
upregulation  of  CRH  release. 

CRH  secretion  activity  is  regulated  by  synaptic  inputs 
received  by  the  PVN  from  multiple  brain  regions  includ¬ 
ing  limbic  structures  such  as  the  hippocampus  and  the 
amygdala,  that  are  activated  during  stress.  It  has  been 
reported  that  for  certain  types  of  stressors,  these  synaptic 
inputs  are  modulated  by  cortisol  independent  of,  or  par¬ 
allel  to,  its  regulatory  function  on  CRH  synthesis  activity 
[34].  On  the  other  hand,  a  series  of  studies  [35-37]  showed 
that  cortisol  did  not  affect  the  basal  spiking  activity  of  the 
PVN.  We  model  the  overall  synaptic  input,  denoted  by 
I(T)  in  Eq.  1,  as  follows 

I(T)=Ibase+Iext(T),  (7) 

where  /base  and  4xt (T)  represent  the  basal  and  stress- 
dependent  synaptic  input  of  the  PVN,  respectively.  As 
the  effect  of  cortisol  on  the  synaptic  input  during  stress 
is  specific  to  the  type  of  stressor  [38-40],  we  assume 
/ext  (4)  to  be  independent  of  O  for  simplicity  and  gener¬ 
ality.  Possible  implications  of  a  cortisol-dependent  input 
function  Iext(T,  O )  on  model  behavior  will  be  discussed  in 
the  Additional  file  1. 

The  secretion  of  CRH  will  also  depend  upon  the 
amount  of  stored  releasable  CRH,  CS(T),  within  the  neu¬ 
ron  and  inside  the  synaptic  vesicles.  Therefore,  Cs  can 
also  be  factored  into  Eq.  1  through  a  source  term  h(Cs) 
which  describes  the  amount  of  CRH  released  per  unit  of 
action  potential  activity  of  CRH  neurons.  Finally,  it  has 
been  hypothesized  that  CRH  enhances  its  own  release 
[23],  especially  when  external  stressors  are  present.  The 
enhancement  of  CRH  release  by  CRH  is  mediated  by 
activation  of  the  membrane-bound  G-protein-coupled 
receptor  CRHR-1  whose  downstream  signaling  path¬ 
ways  operate  on  timescales  from  milliseconds  to  seconds 
[41,  42].  Thus,  self-upregulation  of  CRH  release  can  be 
modeled  by  including  a  positive  and  increasing  function 
gc(C)  in  the  source  term  in  Eq.  1. 

Combining  all  these  factors  involved  in  regulating  the 
secretion  process,  we  can  rewrite  Eq.  1  by  replacing  fc(0) 
with  h(Cs)gc(C)  as  follows 

^  =  pcI(T)h(Cs)gc(Q  -  dcC.  (8) 

di 

In  this  model  (represented  by  Eqs.  6,  8,  2,  5,  and  4), 
cortisol  no  longer  directly  suppresses  CRH  levels,  rather, 


it  decreases  stored  CRH  availability,  Cs,  through  Eq.  6, 
which  in  turn  decreases  the  secretion  rate  of  CRH.  The 
combination  h(Cs)gc(C )  in  Eq.  8  indicates  the  release  rate 
of  stored  CRH  decreases  when  either  Cs  or  C  decrease. 
We  assume  that  synaptic  inputs  into  the  CRH  neurons 
modulate  the  overall  release  process  with  weight  pc . 


Complete  delay-differential  equation  model 

We  are  now  ready  to  incorporate  the  mechanisms 
described  above  into  a  new,  more  comprehensive  math¬ 
ematical  model  of  the  HPA  axis,  which,  in  summary, 
includes 


(i)  A  delayed  response  of  the  adrenal  cortex  to  cortisol 
(Eq.  5). 

(ii)  A  slow  time-scale  negative  feedback  by  cortisol  on 
CRH  synthesis  (through  the  production  term  Coo(O) 
in  Eq.  6). 

(iii)  A  fast-acting  positive  feedback  of  stored  and 
circulating  CRH  on  CRH  release  (through  the  factor 
h(Cs)gc(C)  of  the  production  term  in  Eq.  8). 

Our  complete  mathematical  model  thus  consists  of 
Eqs.  2,  4,  5,  6,  and  8.  We  henceforth  assume  Ja  (OR,  O )  = 
fA(OR)  depends  on  only  the  cortisol-GR  complex  and  use 
Hill-type  functions  for Ja(OR)  and gR (OR)  [13, 14,  16,  17]. 
Our  full  theory  is  characterized  by  the  following  system  of 
delay  differential  equations: 


dCs  Coo(O)  Cs 
dT  “  7c  ’ 

^  =pcI(T)h(Cs)gc(C ) 


u=PAc(^- 

d  T  F  \Ka  +  OR 
dO 

—  =poA(T  -  Td)  -  doO , 
di 


PrI<1 


dcC , 
—  <4  A, 


Ki  +  (OR)2 


dj^R. 


(9) 

(10) 

(ID 

(12) 

(13) 


The  parameters  I<a,r  represent  the  level  of  A  and  R  at 
which  the  negative  or  positive  effect  are  at  their  half  max¬ 
imum  and  1  —  pr  represents  the  basal  production  rate  for 
GR  when  OR  =  0. 

Of  all  the  processes  modeled,  we  will  see  that  the  slow 
negative  feedback  described  in  Eq.  9  will  be  crucial  in 
mediating  transitions  between  stable  states  of  the  sys¬ 
tem.  The  slow  dynamics  will  allow  state  variables  to  cross 
basins  of  attraction  associated  with  each  of  the  stable 
states. 


Nondimensionalization 

To  simplify  the  further  development  and  analysis  of  our 
model,  we  nondimensionalize  Eqs.  9-13  by  rescaling  all 
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variables  and  parameters  in  a  manner  similar  to  that  of 
Walker  et  al.  [13],  as  explicitly  shown  in  the  Additional 
file  1.  We  find 


dcs  Cqq  (o)  cs 
d  t  tc 

d  c 

—  =  q0I(t)h(cs)gc(c )  -  q2c, 
at 

da  c 

d  t  l+p2(or) 

do 

—  =  a(t  -  td)  -  o, 
at 


-P3*> 


dr  {or)2 
d  t  pd  +  (or)2 


+  P5  -  Per, 


(14) 

(15) 

(16) 

(17) 

(18) 


where  cs,c,a,r,o  are  the  dimensionless  versions  of  the 
original  concentrations  Cs,  C,A,R,  O,  respectively.  The 
dimensionless  delay  in  activation  of  cortisol  production  by 
ACTH  is  now  denoted  t&.  All  dimensionless  parameters 
qupu  td,  and  tc  are  combinations  of  the  physical  parame¬ 
ters  and  are  explicitly  given  in  the  Additional  file  1.  The 
functions  Coo(o),  b(cs),  and  gc(c)  are  dimensionless  ver¬ 
sions  of  Coo(O),  h(Cs),  and gc(C),  respectively,  and  will  be 
chosen  phenomenologically  to  be 


Coo(o)  =c0 o  +e  b< 
h(cs )  =1  -  e~kCs, 


gc(c )  =1  - 


Me 

1  +  (qic)n' 


(19) 


The  form  of  Coq  (o)  is  based  on  the  above-mentioned 
exponential  relation  observed  in  adrenalectomized  rats 
[22,  33].  The  parameters  Cqq  and  b  represent  the  minimum 
dimensionless  level  of  stored  CRH  and  the  decay  rate  of 
the  function,  respectively.  The  function  h(cs)  describes 
how  the  rate  of  CRH  release  increases  with  cs.  Since  the 
amount  of  CRH  packaged  in  release  vesicles  is  likely  reg¬ 
ulated,  we  assume  h(cs)  saturates  at  high  cs.  The  choice  of 
a  decreasing  form  for  Coo(o)  implies  that  increasing  cor¬ 
tisol  levels  will  decrease  the  target  level  (or  production 
rate)  of  cs  in  Eq.  14.  The  reduced  production  of  cs  will 
then  lead  to  a  smaller  h(cs)  and  ultimately  to  a  reduced 
release  source  for  c  (Eq.  15).  As  expected,  the  overall  effect 
of  increasing  cortisol  is  a  decrease  in  the  release  rate  of 
CRH.  Finally,  since  the  upregulation  of  CRH  release  by 
circulating  CRH  is  mediated  by  binding  to  CRH  receptor, 
gc(c)  will  be  chosen  to  be  a  Hill-type  function,  with  Hill- 
exponent  n,  similar  in  form  to  the  function  g/?  (OR)  used  in 
Eqs.  13  and  18.  The  parameter  1  —  pc  represents  the  basal 
release  rate  of  CRH  relative  to  the  maximum  release  rate 
and  q±l  represents  the  normalized  CRH  level  at  which  the 
positive  effect  is  at  half-maximum. 


Fast-slow  variable  separation  and  bistability 

Since  we  assume  the  negative  feedback  effect  of  cortisol 
on  synthesis  of  CRH  operates  over  the  longest  character¬ 
istic  timescale  tc  in  the  problem,  the  full  model  must  be 
studied  across  two  separate  timescales,  a  fast  timescale 
t ,  and  a  slow  timescale  r  =  t/tc  =  st.  The  full  model 
(Eqs.  14-18)  can  be  succinctly  written  in  the  form 


dcs 

dt 

dx 

dt 


=  s(Coc(o)  -  cs), 


=  F  (cs,x), 


(20) 

(21) 


where  x  =  (c,  a ,  o,  r)  is  the  vector  of  fast  dynamical  vari¬ 
ables,  and  F(cs,  x)  denotes  the  right-hand-sides  of  Eqs.  15- 
18.  We  refer  to  the  fast  dynamics  described  by  dx/d t  = 
F(cs,x)  as  a  fast  flow .  In  the  s  0  limit,  it  is  also  easy 
to  see  that  to  lowest  order  cs  is  constant  across  the  fast 
timescale  and  is  a  function  of  only  the  slow  variable  r. 

Under  this  timescale  separation,  the  first  component  of 
Eq.  21  (Eq.  15)  can  be  written  as 


dc 

dt 


=  q(cs(r), I)gc(c)  -  q2c , 


(22) 


where  q(cs( r),/)  =  q0Ih{cs{r))  =  q0I{  1  -  e~kc^x))  is 
a  function  of  cs(r)  and  /.  Since  cs  is  a  function  only  of 
the  slow  timescale  r,  q  can  be  viewed  as  a  bifurcation 
parameter  controlling,  over  short  timescales,  the  fast  flow 
described  by  Eq.  22.  Once  c(t )  quickly  reaches  its  non¬ 
oscillating  quasi-equilibrium  value  defined  by  dc/dt  = 
qgc(c)  —  q2c  =  0,  it  can  be  viewed  as  a  parametric  term  in 
Eq.  16  of  the  pituitary-adrenal  (PA)  subsystem. 

Due  to  the  nonlinearity  of  gc(c),  the  equilibrium  value 
c(q)  satisfying  qgc(c)  =  q2c  may  be  multi-valued  depend¬ 
ing  on  q ,  as  shown  in  Fig.  2a  and  b.  For  certain  values 
of  the  free  parameters,  such  as  n,  1  —  /xc,  and  q\,  bista¬ 
bility  can  emerge  through  a  saddle-node  bifurcation  with 
respect  to  the  bifurcation  parameter  q .  Figure  2b  shows 
the  bifurcation  diagram,  i.e.,  the  nullcline  of  c  defined  by 
qgc(c)  =  q2C . 

For  equilibrium  values  of  c  lying  within  a  certain  range, 
the  PA-subsystem  can  exhibit  a  limit  cycle  in  ( a ,  o,  r)  [13] 
that  we  express  as  (a*(6;  c),  o*(6;  c),  r*(0;  c)),  where  0  = 
27it/tp{c)  is  the  phase  along  the  limit  cycle.  The  dynamics 
of  the  PA-subsystem  depicted  in  Fig.  3  indicate  the  range 
of  c  values  that  admit  limit  cycle  behavior  for  {a,  o,  r), 
while  the  fast  c-nullcline  depicted  in  Fig.  2b  restricts  the 
range  of  bistable  c  values.  Thus,  bistable  states  that  also 
support  oscillating  ( a ,  o,  r)  are  possible  only  for  values  of  c 
that  satisfy  both  criteria. 

Since  in  the  s  0  limit,  circulating  CRH  only  feeds 
forward  into  a,  o,  and  r,  a  complete  description  of  all 
the  fast  variables  can  be  constructed  from  just  c  which 
obeys  Eq.  22.  Therefore,  to  visualize  and  approximate  the 
dynamics  of  the  full  five-dimensional  model,  we  only  need 
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Fig.  2  Nonlinear  gc(c)  and  bistability  of  fast  variables,  a  The  stable  states  of  the  decoupled  system  in  Eq.  22  can  be  visualized  as  the  intersection  of 
the  two  functions  qgc(c )  (dashed  curve)  and  <72 c  (gray  line).  For  a  given  Hill-type  function  gc(c),  Eq.  22  can  admit  one  or  two  stable  states  (solid 
circles),  depending  on  function  parameters.  The  unstable  steady  state  is  indicated  by  the  open  circle,  b  Bifurcation  diagram  of  the  decoupled 
system  (Eq.  22)  with  q  as  the  bifurcation  parameter.  Solid  and  dashed  segments  represent  stable  and  unstable  steady  states  of  the  fast  variables, 
respectively.  L  and  U  label  basins  of  attraction  associated  with  the  lower  and  upper  stable  branches  of  the  c-nullcline.  Left  and  right  bifurcation 
points  (qi,  q_)  and  (q r,  cr)  are  indicated.  Fixed  points  of  c  appear  and  disappear  through  saddle  node  bifurcations  as  q  is  varied  between  q l  and  q r 


to  consider  the  2D  projection  onto  the  fast  c  and  slow  cs 
variable.  A  summary  of  the  time-separated  dynamics  of 
the  variables  in  our  model  is  given  in  Fig.  4. 

To  analyze  the  evolution  of  the  slow  variable  cs( t),  we 
write  our  equations  in  terms  of  r  =  st: 


dcs 

dr 

dx 

£  — 


—  (coo(°)  cs)> 

=  F  (cs,x). 


(23) 

(24) 


In  the  £  0  limit,  the  “outer  solution”  F(cs,x)  ~  0 

simply  constrains  the  system  to  be  on  the  fast  c-nullcline 
defined  by  qgc(c)  =  q2c.  The  slow  evolution  of  cs(r)  along 
the  fast  c-nullcline  depends  on  the  value  of  the  fast  vari¬ 
able  o(t)  through  Coo(o).  To  close  the  slow  flow  subsystem 
for  cs(r),  we  fix  c  to  its  equilibrium  value  as  defined  by  the 


fast  subsystem  and  approximate  Cqo(o(c))  in  Eq.  23  by  its 
period-averaged  value 

(Coo(c))  =  r  Coo{0\9;c))^-  =Coo+  [2n  e~h0*(e-’c)f-. 

Jo  Jo 

(25) 

Since  period-averaged  values  of  0*  increases  with  c, 
(coo(c))  is  a  decreasing  function  of  c  under  physiological 
parameter  regimes.  This  period-averaging  approximation 
allows  us  to  relate  the  evolution  of  cs(r)  in  the  slow  sub¬ 
system  directly  to  c.  The  evolution  of  the  slow  subsystem 
is  approximated  by  the  closed  (cs,  c)  system  of  equations 

^  =  <Coc(c)>-cs,  (26) 

dr 

0  =  q0h(cs)I(t)gc(c)  -  q2c.  (27) 


Fig.  3  Dynamics  of  the  oscillating  PA-subsystem  as  a  function  of  fixed  c.  a  Maximum/minimum  and  period-averaged  values  of  ACTH,  a(t),  as  a 
function  of  circulating  CRH.  b  Maximum/minimum  and  period-averaged  values  of  cortisol  o(f).  Within  physiological  CRH  levels,  ACTH,  GR  (not 
shown),  and  cortisol  oscillate.  The  minima,  maxima,  and  period-averaged  cortisol  levels  typically  increase  with  increasing  c.  The  plot  was  generated 
using  dimensionless  variables  c,  a,  and  o  with  parameter  values  specified  in  [1 3]  and  fd  =  1 .44,  corresponding  to  a  delay  of  =  1 5  min 
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slow  variable  fast  variables  x(t) 


cs(x)  c(t)  a(t)  oft)  r(t) 


2D  system  PA  subsystem 

(non-oscillating)  (oscillating) 

Fig.  4  Classification  of  variables.  Variables  of  the  full  five-dimensional 
model  are  grouped  according  to  their  dynamical  behavior.  cs(r)  is  a 
slow  variable,  while  x(t)  =  (c,  a,  o,  r )  are  fast  variables.  Of  these,  ( a ,  o,  r ) 
form  the  typically  oscillatory  PA-subsystem  that  is  recapitulated  by  c. 
In  the  s  =  1  /fc  1  limit,  the  variable  cs(r)  slowly  relaxes  towards  a 
period-averaged  value  (Coo(o(c))).  Therefore,  the  full  model  can  be 
accurately  described  by  its  projection  onto  the  2D  (cs,  c)  phase  space 


with  (Cqq  (c))  evaluated  in  Eq.  25.  By  self-consistently  solv¬ 
ing  Eqs.  26  and  27,  we  can  estimate  trajectories  of  the  full 
model  when  they  are  near  the  c-nullcline  in  the  2D  ( cs ,  c)- 
subsystem.  We  will  verify  this  in  the  following  section. 

Nullcline  structure  and  projected  dynamics 

The  separation  of  timescales  results  in  a  natural  descrip¬ 
tion  of  the  fast  c-nullcline  in  terms  of  the  parameter  q 
(Fig.  2)  and  the  slow  cs -nullcline  (defined  by  the  rela¬ 
tion  cs  =  (coo(c))  relating  cs  to  c)  in  terms  of  c .  How¬ 
ever,  the  c-nullcline  is  plotted  in  the  (#,  c) -plane  while 
the  cs-nullcline  is  defined  in  the  (c,  cs)-plane.  To  plot  the 
nullclines  together,  we  relate  the  equilibrium  value  of  cs, 
(coo(c)),  to  the  q  coordinate  through  the  monotonic  rela- 
tionship  q(cs)  =  qolh({coo(c )»  =  q0I(  1  -  e~k{c°°(c))) 
and  transform  the  cs  variable  into  the  q  parameter  so 


that  both  nullclines  can  be  plotted  together  in  the  (#,  c)- 
plane.  These  transformed  cs -nullclines  will  be  denoted 
“^-nullclines.” 

We  assume  a  fixed  basal  stress  input  7=1  and  plot 
the  ^-nullclines  in  Fig.  5a  for  increasing  values  of  /<,  the 
parameter  governing  the  sensitivity  of  CRH  release  to 
stored  CRH.  From  the  form  h({Coo(c)))  =  (1  —  e~k(^Coo^)) 
both  the  position  and  the  steepness  of  the  ^-nullcline 
in  (#,  c) -space  depend  strongly  on  k.  Figure  5b  shows  a 
fast  c-nullcline  and  a  slow  ^-nullcline  (transformed  cs- 
nullcline)  intersecting  at  both  stable  branches  of  the  fast 
c-nullcline.  Here,  the  flow  field  indicates  that  the  2D  pro¬ 
jected  trajectory  is  governed  by  fast  flow  over  most  of  the 
(i q ,  c) -space. 

How  the  fast  and  slow  nullclines  intersect  controls  the 
long-term  behavior  of  our  model  in  the  small  s  limit. 
In  general,  the  number  of  allowable  nullcline  intersec¬ 
tions  will  depend  on  input  level  I  and  on  parameters 
(qo, . . .  ,p6>b,k,  n,  /xc,£d)- 

Other  parameters  such  as  qo ,  q\ ,  and  /jlc  appear  directly 
in  the  fast  equation  for  c  and  thus  most  strongly  con¬ 
trol  the  fast  c-nullcline.  Figure  6a  shows  that  for  a  basal 
stress  input  of  I  =  1  and  an  intermediate  value  of  /c,  the 
nullclines  cross  at  both  stable  branches  of  the  fast  subsys¬ 
tem.  As  expected,  numerical  simulations  of  our  full  model 
show  the  fast  variables  (a,  o ,  r)  quickly  reaching  their  oscil¬ 
lating  states  defined  by  the  c-nullcline  while  the  slow 
variable#  =  qoIh(cs)  remains  fairly  constant.  Independent 
of  initial  configurations  that  are  not  near  the  c-nullcline  in 
(#,  c) -space,  trajectories  quickly  jump  to  one  of  the  stable 
branches  of  the  c-nullcline  with  little  motion  towards  the 
#-nullcline,  as  indicated  by  §f  in  Fig.  6a. 

Once  near  the  c-nullcline,  say  when  |F(cs,x)|  £, 

the  trajectories  vary  slowly  according  to  Eq.  23.  Here, 


A  B 


Fig.  5  Slow  and  fast  nullclines  and  overall  flow  field,  a  The  nullcline  of  cs  in  the  s  0  limit  is  defined  by  cs  =  (Coq(c)).To  plot  these  slow  nullclines 
together  with  the  fast  c-nullclines,  we  transform  the  variable  cs  and  represent  it  by  q  through  the  relation  q  =  qoh(cs).  These  transformed  nullclines 
then  become  a  function  of  c  and  can  be  plotted  together  with  the  fast  c-nullclines.  For  each  fixed  value  of  c,  o(t;c)  is  computed  by  employing  a 
built-in  DDE  solver  dde23  in  MATLAB.  The  numerical  solution  is  then  used  to  approximate  (Coo(c))  in  Eq.  25  by  Euler's  method.  The  g-nullcline  shifts 
to  the  right  and  gets  steeper  as  k  increases,  b  The  fast  c-nullcline  defined  by  qgc(c )  =  qjc  (black curve)  is  plotted  together  with  the  slow  cs-nullcline 
plotted  in  the  (q,  c)  plane  ("g-nullcline,"  blue  curve).  Here,  two  intersections  arise  corresponding  to  a  high-cortisol  normal  (N)  stable  state  and  a 
low-cortisol  diseased  (D)  stable  state.  The  flow  vector  field  is  predominantly  aligned  with  the  fast  directions  toward  the  c-nullcline 
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Fig.  6  Equilibria  at  the  intersections  of  nullclines.  a  For  intermediate  values  of  k,  there  are  three  intersections,  two  of  them  representing  stable 
equilibria.  Solid  red  lines  are  projections  of  two  trajectories  of  the  full  model,  with  initial  states  indicated  by  red  dots  and  final  stable  states  shown  by 
black  dots.  The  full  trajectories  approach  the  intersections  of  the  g-nullcline  {blue)  and  c-nullcline  {black),  b  For  large  k  there  is  only  one  intersection 
at  the  upper  branch  of  the  c-nullcline.  Two  trajectories  with  initial  states  near  different  branches  of  the  c-nullcline  both  approach  the  unique 
intersection  (black  dot)  on  the  upper  branch.  The  scenario  shown  here  corresponds  to  a  Type  I  nullcline  structure  as  described  in  the  Additional  file  1 


the  slow  variable  cs  relaxes  to  its  steady  state  value  while 
satisfying  the  constraint  F(cs,x)  &  0.  In  (q,  c)—  space, 
the  system  slowly  slides  along  the  c-nullcline  towards  the 
^-nullcline  (the  §s  paths  in  Fig.  6a).  This  latter  phase 
of  the  evolution  continues  until  the  system  reaches  an 
intersection  of  the  two  nullclines,  indicated  by  the  filled 
dot,  at  which  the  reduced  subsystem  in  cs  and  c  reaches 
equilibrium. 

For  certain  values  of  k  and  if  the  fast  variable  c  is  bistable, 
the  two  nullclines  may  intersect  within  each  of  the  two 
stable  branches  of  the  c-nullcline  and  yield  the  two  dis¬ 
tinct  stable  solutions  shown  in  Fig.  6a.  For  large  /<,  the 
two  nullclines  may  only  intersect  on  one  stable  branch  of 
the  c-nullcline  as  shown  in  Fig.  6b.  Trajectories  that  start 
within  the  basin  of  attraction  of  the  lower  stable  branch 
of  the  c-nullcline  (“initial  state  2”  in  Fig.  6b)  will  stay 
on  this  branch  for  a  long  time  before  eventually  sliding 
off  near  the  bifurcation  point  and  jumping  to  the  upper 
stable  branch.  Thus,  the  long-term  behavior  of  the  full 
model  can  be  described  in  terms  of  the  locations  of  the 
intersections  of  nullclines  of  the  reduced  system. 

Results  and  discussion 

The  dual-nullcline  structure  and  existence  of  multiple  sta¬ 
ble  states  discussed  above  results  from  the  separation  of 
slow  CRH  synthesis  process  and  fast  CRH  secretion  pro¬ 
cess.  This  natural  physiological  separation  of  time  scales 
ultimately  gives  rise  to  slow  dynamics  along  the  fast  c- 
nullcline  during  stress.  The  extent  of  this  slow  dynamics 
will  ultimately  determine  whether  a  transition  between 
stable  states  can  be  induced  by  stress.  In  this  section, 
we  explore  how  external  stress-driven  transitions  medi¬ 
ated  by  the  fast-slow  negative  feedback  depend  on  system 
parameters. 


Changes  in  parameters  that  accompany  trauma  can  lead 
to  shifts  in  the  position  of  the  nullclines.  For  example, 
if  the  stored  CRH  release  process  is  sufficiently  compro¬ 
mised  by  trauma  (smaller  /<),  the  slow  ^-nullcline  moves  to 
the  left,  driving  a  bistable  or  fully  resistant  organism  into  a 
stable  diseased  state.  Interventions  that  increase  k  would 
need  to  overcome  hysteresis  in  order  to  restore  normal 
HPA  function.  More  permanent  changes  in  parameters 
are  likely  to  be  caused  by  physical  rather  than  by  psycho¬ 
logical  traumas  since  such  changes  would  imply  altered 
physiology  and  biochemistry  of  the  person.  Traumatic 
brain  injury  (TBI)  is  an  example  in  which  parameters  can 
be  changed  permanently  by  physical  trauma.  The  injury 
may  decrease  the  sensitivity  of  the  pituitary  to  cortisol- 
GR  complex,  which  can  be  described  by  decreasing  p2  in 
our  model.  Such  change  in  parameter  would  lead  to  a  left¬ 
ward  shift  of  the  ^-nullcline  and  an  increased  likelihood  of 
hypocortisolism. 

In  the  remainder  of  this  work,  we  focus  on  how  exter¬ 
nal  stress  inputs  can  by  themselves  induce  stable  but 
reversible  transitions  in  HPA  dynamics  without  changes 
in  physiological  parameters.  Specifically,  we  consider 
only  temporary  changes  in  I(t)  and  consider  the  time- 
autonomous  problem.  Since  the  majority  of  neural  circuits 
that  project  to  the  PVN  are  excitatory  [43],  we  assume 
external  stress  stimulates  CRH  neurons  to  release  CRH 
above  its  unit  basal  rate  and  that  I{t)  =  1  +  lext(t) 
(4ase  =  1)  With  /ext  >  0. 

To  be  more  concrete  in  our  analysis,  we  now  choose 
our  nullclines  by  specifying  parameter  values.  We  esti¬ 
mate  many  of  the  dimensionless  parameters  by  using 
values  from  previous  studies,  as  listed  in  Table  1.  Of  the 
four  remaining  parameters,  /xc,  qo ,  q\,  and  /<,  we  will  study 
how  our  model  depends  on  k  while  fixing  /Ac,qo,  and 
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Table  1  Dimensionless  parameter  values  of  our  full  model 


Parameter 

Value 

Source  and  Ref. 

Description 

n 

5 

Assumed 

Hill  coefficient  in  upregulation  function  gc(c) 

Coo 

0.2 

Estimated  from  [22] 

Baseline  stored  CRH  level 

b 

0.6 

Estimated  from  [22] 

Relates  cortisol  to  stored  CRH  level 

k 

Undetermined 

Relates  stored  CRH  to  CRH  release  rate 

ilc 

Undetermined 

Basal  CRH  release  rate 

qo 

Undetermined 

Maximum  CRH  release  rate 

Undetermined 

Circulating  CRH  for  half-maximum  self-upregulation 

Q2 

1.8 

Estimated  from  [21] 

Ratio  of  CRH  and  cortisol  decay  rates 

Pl] 

0.067 

P2]  [H] 

(or)-complex  level  for  half-maximum  feedback 

Ps 

7.2 

Ps  [13] 

Ratio  of  ACTH  and  cortisol  decay  rates 

Pa 

0.05 

Pa  [13] 

(or)-complex  level  for  half-maximum  upregulation 

Ps 

0.11 

Ps  [13] 

Basal  GR  production  rate  by  pituitary 

P6 

2.9 

Pe  [13] 

Ratio  of  GR  and  cortisol  decay  rates 

k 

69.3 

Assumed 

CRH  biosynthesis  timescale 

k 

1.44 

V[13] 

Delay  in  ACTH-activated  cortisol  release 

q\.  Three  possible  nullcline  configurations  arise  accord- 
ing  to  the  values  of  ptC)  qo,  and  q\  and  are  delineated  in 
the  Additional  file  1.  We  have  also  implicitly  considered 
only  parameter  regimes  that  yield  oscillations  in  the  PA 
subsystem  at  the  stable  states  defined  by  the  nullcline 
intersections. 

Given  these  considerations,  we  henceforth  chose  ptc  = 
0.6,  q\  =  0.04,  and  qo  =  77.8  for  the  rest  of  our  analysis. 
This  choice  of  parameters  is  motivated  in  the  Additional 
file  1  and  corresponds  to  a  so-called  “Type  I”  nullcline 
structure.  In  this  case,  three  possibilities  arise:  one  inter¬ 
section  on  the  lower  stable  branch  of  the  c-nullcline  if 
k  <  /cl,  two  intersections  if  /<l  <  k  <  /<r  (Fig.  6a),  and  one 
intersection  on  the  upper  stable  branch  of  the  c-nullcline 


if  k  >  /<r  (Fig.  6b).  For  our  chosen  set  of  parameters  and 
a  basal  stress  input  7=1,  the  critical  values  /<l  =  2.50  < 
/cr  =  2.54  are  given  by  Eq.  A3  in  the  Additional  file  1. 

Normal  stress  response 

Activation  of  the  HPA  axis  by  acute  stress  culminates  in 
an  increased  secretion  of  all  three  main  hormones  of  the 
HPA  axis.  Persistent  hypersecretion  may  lead  to  numer¬ 
ous  metabolic,  affective,  and  psychotic  dysfunctions 
[44,  45].  Therefore,  recovery  after  stress-induced  pertur¬ 
bation  is  essential  to  normal  HPA  function.  We  explore 
the  stability  of  the  HPA  axis  by  initiating  the  system  in  the 
upper  of  the  two  stable  points  shown  in  Fig.  7a  and  impos¬ 
ing  a  120  min  external  stress  input  of/ext  =  0.1.  The  HPA 


A  B 


Fig.  7  Normal  stress  response.  Numerical  solution  for  the  response  to  a  1 20  min  external  stress  /ex t  =  0.1 .  a  At  the  moment  the  external  stress  is 
turned  on,  the  value  of  (q,c)  increases  from  its  initial  stable  solution  at  (64.4,27)  to  (71,27)  after  which  the  circulating  CRH  level  c,  quickly  reaches 
the  fast  c-nullcline  (block)  before  slowly  evolving  along  it  towards  the  slow  g-nullcline  (blue).  After  short  durations  of  stress,  the  system  returns  to  its 
starting  point  within  the  normal  state  basin,  b  The  peaks  of  the  cortisol  level  are  increased  during  stress  (red)  but  return  to  their  original  oscillating 
values  after  the  stress  is  turned  off 
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axis  responds  with  an  increase  in  the  peak  level  of  cortisol 
before  relaxing  back  to  its  original  state  after  stress  is  ter¬ 
minated  (Fig.  7b).  This  transient  process  is  depicted  in  the 
projected  ( q ,  c) -space  in  Fig.  7a. 

Upon  turning  on  stress,  the  lumped  parameter  q  and 
the  slow  nullcline  shift  to  the  right  by  10  %  since  q  = 
#o(l+4xt)^((coo(c)))  (see  Fig.  7a).  The  trajectory  will  then 
move  rapidly  upward  towards  the  new  value  of  c  on  the 
c-nullcline;  afterwards,  it  moves  very  slowly  along  the  c- 
nullcline  towards  the  shifted  nullcline.  After  120  min, 
the  system  arrives  at  the  “x”  on  the  c-nullcline  (Fig.  7a). 
Once  the  stress  is  shut  off  the  ^-nullcline  returns  to  its 
original  position  defined  by  /  =  1.  The  trajectory  also 
jumps  horizontally  back  to  near  the  initial  q  value  and 
quickly  returns  to  the  original  upper-branch  stable  point. 

External  stress  induces  transition  from  normal  to  diseased 
state 

We  now  discuss  how  transitions  from  a  normal  to  a  dis¬ 
eased  state  can  be  induced  by  positive  (excitatory)  external 
stress  of  sufficient  duration.  In  Fig.  8,  we  start  the  system 
in  the  normal  high-c  state. 

Upon  stimulation  of  the  CRH  neurons  through  7ext  >  0, 
both  CRH  and  average  glucocorticoid  levels  are  increased 


while  the  average  value  of  Coo(o(t))  is  decreased  since 
Coo(o)  is  a  decreasing  function  of  o .  As  cs(r)  slowly  decays 
towards  the  decreased  target  value  of  (Coo(o(c))),  h(cs(z)) 
and  hence  q(cs),  also  decrease.  As  shown  in  Fig.  8a,  much 
of  this  decrease  occurs  along  the  high-c  stable  branch  of 
the  c-nullcline.  Once  the  external  stress  is  switched  off, 
q  will  jump  back  down  by  a  factor  of  1/(1  +  7ext).  If  the 
net  decrease  in  q  is  sufficient  to  bring  it  below  the  bifur¬ 
cation  value  qi  &  64  at  the  leftmost  point  of  the  upper 
knee,  the  system  crosses  the  separatrix  and  approaches 
the  alternate,  diseased  state.  Thus,  the  normal-to-diseased 
transition  is  more  likely  to  occur  if  the  external  stress  is 
maintained  long  enough  to  cause  a  large  net  decrease  in 
q,  which  includes  the  decrease  in  q  incurred  during  the 
slow  relaxation  phase,  plus  the  drop  in  q  associated  with 
cessation  of  stress.  The  minimum  duration  required  for 
normal-to-diseased  transition  should  also  depend  on  the 
magnitude  of  Iext.  The  relation  between  the  stressor  mag¬ 
nitude  and  duration  will  be  illustrated  in  the  Additional 
file  1  (Fig.  A4). 

A  numerical  simulation  of  our  model  with  a  30  hr  Iext  = 
0.2  was  performed,  and  the  trajectory  in  (q,  c)-space  is 
shown  in  Fig.  8a.  The  corresponding  cortisol  level  along 
this  trajectory  is  plotted  in  Fig.  8b,  showing  that  indeed 


Fig.  8  Stress-induced  transitions  into  an  oscillating  low-cortisol  diseased  state.  An  excitatory  external  stress  /ex t  =  0.2  is  applied  for  30  hrs.  Here,  the 
system  reaches  the  new  stable  point  set  by  /  =  1.2  before  stress  is  terminated  and  the  g-nullcline  reverts  to  its  original  position  set  by  /  =  1 .  a  At 
intermediate  values  of  2.50  <  k  <  2.54,  when  two  stable  state  arise,  a  transition  from  the  normal  high-cortisol  state  into  the  diseased  low-cortisol 
state  can  be  induced  by  chronic  external  stress,  b  Numerical  solutions  of  cortisol  level  o(T )  plotted  against  the  original  time  variable  T  shows  the 
transition  to  the  low-cortisol  diseased  state  shortly  after  cessation  of  stress,  c  and  d  If  k  >  /cr  =  2.54,  only  the  normal  stable  state  exists.  The  system 
will  recover  and  return  to  its  original  normal  state  after  a  transient  period  of  low  cortisol 
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a  stable  transition  to  the  lower  cortisol  state  occurred 
shortly  after  the  cessation  of  stress. 

In  addition  to  a  long-term  external  stress,  the  stable 
transition  to  a  diseased  state  requires  2.50  <  k  <  2.54 
and  the  existence  of  two  stable  points.  On  the  other  hand, 
when  k  >  /cr  =  2.54,  the  enhanced  CRH  release  stim¬ 
ulates  enough  cortisol  production  to  drive  the  sole  long 
term  solution  to  the  stable  upper  normal  branch  of  the 
c-nullcline,  rendering  the  HPA  system  resistant  to  stress- 
induced  transitions. 

The  response  to  chronic  stress  initially  follows  the  same 
pattern  as  described  above  for  the  two -stable -state  case, 
as  shown  in  Fig.  8c.  However,  the  system  will  continue 
to  evolve  along  the  lower  branch  towards  the  #-nullcline, 
eventually  sliding  off  the  lower  branch  near  the  right 
bifurcation  point  (indicated  in  Fig.  2b  and  Fig.  A2  in  the 
Additional  file  1  by  (#r,  cr))  and  returning  to  the  sin¬ 
gle  normal  equilibrium  state.  Thus,  when  k  is  sufficiently 
high,  the  system  may  experience  a  transient  period  of 
lowered  cortisol  levels  after  chronic  stress  but  will  even¬ 
tually  recover  and  return  to  the  normal  cortisol  state.  The 
corresponding  cortisol  level  shown  in  Fig.  8d  shows  this 
recovery  at  T  &  3400  min,  which  occurs  approximately 
1500  min  after  the  cessation  of  stress. 


Transition  to  diseased  state  depends  on  stress  timing 

We  have  shown  how  transitions  between  the  oscillating 
normal  and  diseased  states  depend  on  the  duration  of  the 
external  stress  7ext-  However,  whether  a  transition  occurs 
also  depends  on  the  time  -  relative  to  the  phase  of  the 
intrinsic  ultradian  oscillations  -  at  which  a  fixed- duration 
external  stress  is  initiated.  To  illustrate  this  dependence 
on  phase,  we  plot  in  Fig.  9a  and  b  two  solutions  for  o(T ) 
obtained  with  a  250  min  Iext  =  0.1  initiated  at  differ¬ 
ent  phases  of  the  underlying  cortisol  oscillation.  If  stress 
is  initiated  near  the  nadir  of  the  oscillations,  a  transition 
to  the  low-cortisol  diseased  state  occurs  and  is  com¬ 
pleted  at  approximately  T  =  1000  min  (Fig.  9a,  c).  If, 
however,  stress  is  initiated  near  the  peak  of  the  oscil¬ 
lations,  the  transition  does  not  occur  and  the  system 
returns  to  the  normal  stable  state  (Fig.  9b,  d).  In  this 
case,  a  longer  stress  duration  would  be  required  to  push 
the  trajectory  past  the  low-#  separatrix  into  the  diseased 
state. 

As  discussed  earlier,  an  increase  in  period-averaged  cor¬ 
tisol  level  during  stress  drives  a  normal-to-diseased  state 
transition.  We  see  that  the  period-averaged  level  of  corti¬ 
sol  under  increased  stress  is  different  for  stress  started  at 
120  min  from  stress  started  at  150  min. 


A 


B 


b 


c 


Fig.  9  Stress  tinning  and  transition  to  low-cortisol  oscillating  state.  Cortisol  levels  in  response  to  /ex t  =  0.1  applied  over  250  min.  a  If  stress  is  initiated 
at  T  =  1 50  min,  a  transition  to  the  low-cortisol  diseased  state  is  triggered,  b  If  stress  is  initiated  at  T  =  1 20  min,  the  system  returns  to  its  normal 
high-cortisol  state.  Note  that  the  first  peak  (marked  by  "▼")  during  the  stress  in  (a)  is  higher  than  the  first  peak  in  (b).  c  If  stress  is  initiated  at 
T  =  150  min,  stress  cessation  and  the  slow  relaxation  along  the  c-nullcline  during  stress  are  sufficient  to  bring  gjust  left  of  the  separatrix,  inducing 
the  transition,  d  For  initiation  time  T  =  120  min,  q  remains  to  the  right  of  the  separatrix,  precluding  the  transition 


Ki  m  et  al.  Biology  Direct  (20 1 6)  1 1 : 1 3 


Page  13  of  18 


As  detailed  in  the  Additional  file  1  (Fig.  A5),  the  ampli¬ 
tude  of  the  first  cortisol  peak  after  the  start  of  stress  is 
significantly  lower  when  the  applied  stress  is  started  dur¬ 
ing  the  falling  phase  of  the  intrinsic  cortisol  oscillations. 
The  difference  between  initial  responses  in  o(t )  affects 
the  period-averaging  in  (Coo(o))  during  external  stress, 
ultimately  influencing  cs  and  consequently  determining 
whether  a  transition  occurs.  Note  that  this  phase  depen¬ 
dence  is  appreciable  only  when  stress  duration  is  near  the 
threshold  value  that  brings  the  system  close  to  the  sepa- 
ratrix  between  normal  and  diseased  basins  of  attraction. 
Trajectories  that  pass  near  separatrices  are  sensitive  to 
small  changes  in  the  overall  negative  feedback  of  cortisol 
on  CRH  synthesis,  which  depend  on  the  start  time  of  the 
stress  signal. 

Stress  of  intermediate  duration  can  induce  "reverse" 
transitions 

We  can  now  use  our  theory  to  study  how  positive  stressors 
/ext  may  be  used  to  induce  “reverse”  transitions  from  the 
diseased  to  the  normal  state.  Understanding  these  reverse 
transitions  may  be  very  useful  in  the  context  of  exposure 
therapy  (ET),  where  PTSD  patients  are  subjected  to  stres¬ 
sors  in  a  controlled  and  safe  manner,  using  for  example, 
computer-simulated  “virtual  reality  exposure.”  Within  our 


model  we  can  describe  ET  as  external  stress  (/ex t  >  0) 
applied  to  a  system  in  the  stable  low-c  diseased  state.  The 
resulting  horizontal  shift  in  q  causes  the  system  to  move 
rightward  across  the  separatrix  and  suggests  a  transition 
to  the  high-c  normal  state  can  occur  upon  termination  of 
stress. 

As  shown  in  Fig.  10a,  if  stressor  of  sufficient  duration 
is  applied,  the  trajectory  reaches  a  point  above  the  unsta¬ 
ble  branch  of  the  c-nullcline  upon  termination  leading  to 
the  normal,  high-cortisol  state  (Fig.  10b).  Since  the  initial 
motion  is  governed  by  fast  flow,  the  minimum  stress  dura¬ 
tion  needed  to  incite  the  diseased-to-normal  transition  is 
short,  on  the  timescale  of  minutes.  However,  if  the  stressor 
is  applied  for  too  long,  a  large  reduction  in  q  is  experi¬ 
enced  along  the  upper  stable  branch.  Cessation  of  stress 
might  then  lower  q  back  into  the  basin  of  attraction  of 
the  low-cortisol  diseased  state  (Fig.  10c).  Figure  lOd  shows 
the  cortisol  level  transiently  increasing  to  a  normal  level 
before  reverting  back  to  low  levels  after  approximately 
1400  min. 

Within  our  dynamical  model,  stresses  need  to  be  of 
intermediate  duration  in  order  to  induce  a  stable  transi¬ 
tion  from  the  diseased  to  the  normal  state.  The  occurrence 
of  a  reverse  transition  may  also  depend  on  the  phase  (rel¬ 
ative  to  the  intrinsic  oscillations  of  the  fast  PA  subsystem) 
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Fig.  10  Stress-induced  transitions  to  high-cortisol  oscillating  state,  a  Projected  2D  system  dynamics  when  a  stressor  of  amplitude  /ex t  =  0.2  is 
applied  for  9  min  starting  at  T  =  1 20  min.  c  is  increased  just  above  the  unstable  branch  (c  «  20)  to  allow  the  unstressed  system  to  cross  the 
separatrix  and  transition  to  the  normal  high-c  stable  state.  bThe  plot  of  o(7~)  shows  the  transition  to  the  high-cortisol,  high-oscillation  amplitude 
state  shortly  after  the  9  min  stress,  c  A  stressor  turned  off  after  780  min  (13  hrs)  leaves  the  system  in  the  basin  of  attraction  of  the  diseased  state, 
d  Cortisol  levels  are  pushed  up  but  after  about  1400  min  relax  back  to  levels  of  the  original  diseased  state 
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over  which  stress  was  applied,  especially  when  the  stress 
duration  is  near  its  transition  thresholds.  For  a  reverse 
diseased-to-normal  transition  to  occur,  the  decrease  in  cs 
cannot  be  so  large  that  it  brings  the  trajectory  past  the  left 
separatrix,  as  shown  in  Fig.  10c.  Therefore,  near  the  max¬ 
imum  duration,  stress  initiated  near  the  nadir  of  cortisol 
oscillation  will  be  more  effective  at  triggering  the  transi¬ 
tion  to  a  normal  high-cortisol  state.  Overall,  these  results 
imply  that  exposure  therapy  may  be  tuned  to  drive  the 
dynamics  of  the  HPA  axis  to  a  normal  state  in  patients 
with  hypocortisolism-associated  stress  disorders. 

Conclusions 

We  developed  a  theory  of  HPA  dynamics  that  includes 
stored  CRH,  circulating  CRH,  ACTH,  cortisol  and  glu¬ 
cocorticoid  receptor.  Our  model  incorporates  a  fast  self- 
upregulation  of  CRH  release,  a  slow  negative  feedback 
effect  of  cortisol  on  CRH  synthesis,  and  a  delay  in  ACTH- 
activated  cortisol  synthesis.  These  ingredients  allow  our 
model  to  be  separated  into  slow  and  fast  components  and 
projected  on  a  2D  subspace  for  analysis. 

Depending  on  physiological  parameter  values,  there 
may  exist  zero,  one,  or  two  stable  simultaneous  solutions 
of  both  fast  and  slow  variables.  For  small  /<,  the  param¬ 
eter  that  relates  the  amount  of  stored  releasable  CRH 
vesicles  to  its  secretion  rate,  CRH  release  is  weak  and 
only  the  low-CRH  equilibrium  point  arises;  an  individ¬ 
ual  with  such  k  is  trapped  in  the  low-cortisol  “diseased” 
state.  For  large  /<,  only  the  high-CRH  normal  state  arises, 
rendering  the  individual  resistant  to  acquiring  the  long¬ 
term,  low-cortisol  side-effect  of  certain  stress  disorders. 
When  only  one  stable  solution  arises,  HPA  dysregulation 
must  depend  on  changes  in  parameters  resulting  from 
permanent  physiological  modifications  due  to  e.g.,  aging, 
physical  trauma,  or  stress  itself  [45,  46].  For  example,  it 
has  been  observed  that  older  rats  exhibit  increased  CRH 
secretion  while  maintaining  normal  levels  of  CRH  mRNA 
in  the  PVN  [47].  Such  a  change  could  be  interpreted  as  an 
age-dependent  increase  in  /c,  which,  in  our  model,  implies 
that  aging  makes  the  organism  more  resistant  to  stress- 
induced  hypocortisolism.  Indeed,  it  has  been  suggested 
that  prevalence  of  PTSD  declines  with  age  [48,  49]. 

Other  regulatory  systems  that  interact  with  or  regu¬ 
late  the  HPA  axis  can  also  affect  parameter  values  in 
our  model.  Gonadal  steroids,  which  are  regulated  by 
another  neuroendocrine  system  called  the  hypothalamic- 
pituitary-gonadal  (HPG)  axis,  activate  the  preoptic  area 
(POA)  of  the  hypothalamus  [50,  51],  which  in  turn  atten¬ 
uates  the  excitatory  effects  of  medial  amygdala  stimula¬ 
tion  of  the  HPA  axis  [52].  Thus,  low  testosterone  levels 
associated  with  hypogonadism  would  effectively  increase 
I {t)  within  our  model,  shift  the  ^-nullcline  in  the  (q,  c)~ 
space,  and  in  turn  increase  cortisol  levels.  One  might 
consider  this  as  a  possible  explanation  for  chronically 


elevated  cortisol  levels  observed  in  major  depressive  dis¬ 
order  patients  who  suffers  from  hypogonadism.  Although 
it  is  beyond  the  scope  of  this  paper,  one  may  further 
investigate  role  of  gonadal  hormones,  or  role  of  any 
other  interacting  systems,  in  mediating  stress  response  by 
considering  which  parameters  would  be  affected  in  our 
model. 

Within  certain  parameter  regimes  and  for  intermedi¬ 
ate  k ,  our  theory  can  also  exhibit  bistability.  When  two 
stable  solutions  arise,  we  identify  the  states  with  low  oscil¬ 
lating  levels  of  cortisol  as  the  diseased  state  associated 
with  hypocortisolism.  Transitions  between  different  sta¬ 
ble  states  can  be  induced  by  temporary  external  stress 
inputs,  implying  that  HPA  dysregulation  may  develop 
without  permanent  “structural”  or  physiological  changes. 
Stresses  that  affect  secretion  of  CRH  by  the  PVN  are 
shown  to  be  capable  of  inducing  transitions  from  normal 
to  diseased  states  provided  they  are  of  sufficient  duration 
(Fig.  8). 

Our  model  offers  a  mechanistic  explanation  to  the 
seemingly  counter-intuitive  phenomenon  of  lower  corti¬ 
sol  levels  after  stress-induced  activation  of  cortisol  pro¬ 
duction.  Solutions  to  our  model  demonstrate  that  the 
negative-feedback  effect  of  a  temporary  increase  in  corti¬ 
sol  on  the  synthesis  process  of  CRH  can  slowly  accumulate 
during  the  stress  response  and  eventually  shift  the  system 
into  a  different  basin  of  attraction.  Such  a  mechanism  pro¬ 
vides  an  alternative  to  the  hypothesis  that  hypocortisolism 
in  PTSD  patients  results  from  permanent  changes  in  phys¬ 
iological  parameters  associated  with  negative-feedback  of 
cortisol  [53,  54]. 

We  also  find  that  external  stress  can  induce  the  “reverse” 
transition  from  a  diseased  low-cortisol  state  to  the  nor¬ 
mal  high-cortisol  state.  Our  results  imply  that  re-exposure 
to  stresses  of  intermediate  duration  can  drive  the  system 
back  to  normal  HPA  function,  possibly  “decoupling”  stress 
disorders  from  hypocortisolism. 

Interestingly,  we  show  that  the  minimum  duration 
required  for  either  type  of  transition  depends  on  the  time 
at  which  the  stress  is  initiated  relative  to  the  phase  of  the 
intrinsic  oscillations  in  (a,o,r).  Due  to  subtle  differences 
in  cortisol  levels  immediately  following  stress  initiation 
at  different  phases  of  the  intrinsic  cortisol  oscillation,  the 
different  cumulative  negative-feedback  effect  on  CRH  can 
determine  whether  or  not  a  trajectory  crosses  the  separa¬ 
trix  (Fig.  9).  When  the  duration  of  external  stress  is  near  its 
threshold,  normal-to-diseased  state  transitions  are  easier 
to  induce  when  stress  is  initiated  near  the  nadir  of  cor¬ 
tisol  oscillations.  Reverse  diseased-to-normal  transitions 
are  more  easily  induced  when  stress  is  initiated  near  the 
peak  of  the  oscillations. 

In  summary,  our  theory  provides  a  mechanistic  picture 
that  connects  ortisol  dysregulation  with  stress  disorders 
and  a  mathematical  framework  one  can  use  to  study  the 
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downstream  effects  of  therapies  such  as  brief  eclectic  psy¬ 
chotherapy  (BEP)  and  exposure  therapy  (ET).  Both  thera¬ 
pies  involve  re-experiencing  stressful  situations  directly  or 
through  imagination,  and  have  been  consistently  proven 
effective  as  first-line  treatments  for  PTSD  symptoms 
[55-57].  Our  results  suggest  that  ET  can  directly  alter  and 
“decouple”  the  expression  of  cortisol  from  an  underlying 
upstream  disorder.  Changes  in  neuronal  wiring  that  typi¬ 
cally  occur  over  slower  times  scales  is  also  expected  after 
ET  [58].  In  our  model,  such  changes  would  lead  to  slow 
variations  in  the  basal  input  1(f).  Thus,  cortisol  level  may 
not  be  tightly  correlated  with  PTSD,  particularly  in  the 
context  of  ET. 

It  is  important  to  emphasize  that  we  modeled  neu¬ 
roendocrine  dynamics  downstream  of  the  stress  input 
7ext-  Understanding  how  the  form  of  the  stress  function 
/ext  depends  on  the  type  of  stress  experienced  requires  a 
more  detailed  study  of  more  upstream  processes,  includ¬ 
ing  how  hormones  might  feed  back  on  these  higher-brain 
processes.  Since  higher  cortisol  levels  are  found  among 
female  PTSD  patients  with  a  history  of  childhood  abuse 
[59]  and  among  PTSD  patients  who  have  experienced  a 
nuclear  accident  [60],  future  studies  of  such  divergent, 
experience-dependent  dysregulation  will  rely  on  more 
complex  input  functions  /ext(0-  For  example,  under  peri¬ 
odic  driving,  complex  resonant  behavior  should  arise 
depending  on  the  amplitude  and  frequency  of  the  exter¬ 
nal  stress  /ex tif)  and  the  nullcline  structure  of  the  specific 
system.  Moreover,  effects  of  other  regulatory  networks 
that  interacts  with  the  HPA  axis  can  be  included  in  our 
model  through  appropriate  forms  of  /ext(0-  For  exam¬ 
ple,  the  effects  of  gonadal  steroids  in  the  stress  response 
mentioned  above  can  be  further  investigated  by  con¬ 
sidering  a  form  of  /ext(0  that  is  dependent  on  gonadal 
steroids  level.  Many  other  interesting  properties,  such  as 
response  to  dexamethasone  administration,  can  be  read¬ 
ily  investigated  within  our  model  under  different  system 
parameters. 

Reviewers'  comments 

Reviewer's  report  1 :  Daniel  Coombs,  University  of  British 
Columbia,  Canada 

Summary:  I  find  this  to  be  a  very  interesting  extension 
of  previous  modeling  efforts  on  an  important  neuroen¬ 
docrine  system.  This  is  a  great  example  of  dynamical 
systems  application  to  physiology  and  the  interpretation 
of  normal  and  disease  states  as  attractors  of  the  system. 
The  findings  are  based  squarely  on  the  modeling  and 
the  authors  related  their  model  and  findings  to  experi¬ 
mental  data  and  propose  extensions  to  include  additional 
biological  knowledge  in  the  model.  Potential  caveats  are 
uncertainties  in  parameter  estimates,  omission  of  impor¬ 
tant  interactions  with  other  physiological  systems,  and 


some  question  of  how  external  forcing  (positive  and  neg¬ 
ative  stressors)  should  be  input  to  the  dynamical  system. 
These  possible  shortcomings  are  acknowledged  in  the 
manuscript  and  provide  motivation  for  future  work  within 
the  existing  framework. 

Authors’  response:  We  thank  Dr  Daniel  Coombs  for  his 
appreciation  of  the  work. 

Recommendations:  Many  of  the  results  appear  to 
depend  sensitively  on  the  parameter  k.  This  parameter 
describes  the  regulation  of  CRH  secretion  by  stored  CRH 
(mathematically,  the  link  between  dC/dt  and  Cs).  For 
example,  on  page  16,  it  is  stated  that  you  need  2.5  <  k  < 
2.54  to  find  a  certain  transition  to  the  disease  state  of  the 
model.  The  need  for  precision  of  a  few  %  in  a  param¬ 
eter  value  might  suggest  that  the  system  is  not  terribly 
robust.  Can  you  comment  on  this?  Is  there  any  basis  on 
which  we  might  estimate  k  for  a  given  individual?  Could 
k  be  manipulated  somehow?  Also,  for  readability:  I  think 
it  would  be  good  to  remind  the  reader  of  the  precise 
biological  meaning  of  this  important  parameter  before 
describing  its  effects  in  the  model  (in  the  “Conclusions” 
section). 

Authors’  response:  We  thank  Dr.  Coombs  for  raising  the 
important  issue  of  robustness  of  dynamical  structure  of  our 
model.  The  estimated  range  of  0.04  ofk  (2.5  <  k  <  2.54) for 
bistability  may  be  interpreted  too  narrow  or  wide  depend¬ 
ing  on  the  range  of  values  ofk  observed  from  experimental 
data.  To  best  of  our  knowledge ,  such  data  does  not  yet 
exist.  Although  direct  measurements  on  humans  are  diffi¬ 
cult ,  there  have  been  in  vitro  studies  on  rats  that  estimated 
the  releasable  pool  size  of  synaptic  vesicles  from  hippocam¬ 
pal  neurons  [61].  We  believe  one  could  design  a  similar 
experiment  to  estimate  the  size  of  stored  releasable  CRH  in 
the  PVN.  We  are  also  very  much  interested  in  knowing  ifk 
can  be  modulated  and  manipulated ,  as  it  will  have  impor¬ 
tant  implications  on  our  model  prediction  of  disease  onset 
and  development. 

The  range  of  k  used  for  analysis  in  the  manuscript 
was  chosen  arbitrarily  for  the  purpose  of  concreteness. 
Depending  on  other  parameters  qo ,  q\,  and  pc,  the  range 
ofk  could  be  significantly  greater.  In  particular,  the  range 
of  possible  k  that  give  rise  to  bistability  is  infinite  in 
parameter  regimes  that  correspond  to  “ Type  IT'  null¬ 
cline  structure.  We  have  chosen  a  k  value  that  corre¬ 
spond  to  a  “ Type  I"  nullcline  for  for  the  richness  of  the 
corresponding  mathematical  structure.  Details  on  criti¬ 
cal  values  1<l,r  and  possible  nullcline  configurations  are 
discussed  and  illustrated  (Fig.  A2)  in  the  Additional 
file  1. 

We  also  thank  Dr.  Coombs  for  his  suggestion  to  recall  the 
biological  meaning  of  the  parameter  in  the  Summary  and 
Conclusion  section.  We  have  edited  the  section  accordingly. 

Minor  issues:  No  minor  issues,  the  paper  is  written  well 
and  the  figures  are  clear. 
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Reviewer's  report  2:  Yang  Kuang,  Arizona  State  University, 
United  States  of  America 

Summary:  The  authors  carefully  extended  some  existing 
models  of  oscillating  neuroendocrine  dynamics  by  explic¬ 
itly  incorporating,  in  the  cortisol  equation,  a  discrete 
time  delay  Ta  representing  the  time  needed  for  ACTH- 
mediated  activation  of  cortisol  production  downstream 
of  the  hypothalamus.  Based  on  a  systematical  analysis 
of  this  more  plausible  model,  they  developed  a  com¬ 
pelling  theory  of  the  HPA  dynamics  that  includes  stored 
CRH,  circulating  CRH,  ACTH,  cortisol  and  glucocorti¬ 
coid  receptor.  This  is  a  significant  enhancement  over  the 
existing  theories. 

Authors'  response:  We  thank  Dr.  Yang  Kuang  for  a  posi¬ 
tive  summary  of  our  manuscript. 

Recommendation:  To  facilitate  a  deeper  biological 
appreciation  of  the  fast  and  slow  dynamics  described  by 
the  model,  I  suggest  the  authors  also  present  a  parameter 
table  for  the  original  parameters,  including  their  values, 
units  and  references.  Some  details  on  estimation  of  these 
parameters  are  also  helpful. 

Authors'  response:  We  thank  Dr.  Kuang  for  his  sugges¬ 
tions.  The  list  of  parameters  used  in  our  analysis  and 
numerical  simulations  are  now  listed  in  Table  1  in  the 
manuscript ,  with  details  on  their  estimations  included  in 
the  Additional  file  1. 

Minor  issues:  None. 

Reviewer's  report  3:  Ha  Youn  Lee,  Keck  School  of  Medicine 
University  of  Southern  California,  United  States  of  America 
Summary:  Kim  et  al.  modified  the  existing  cortisol 
dynamics  model  in  the  HPA  axis,  mainly  to  replicate 
stable  ultradian  oscillations  in  cortisol.  By  introducing 
the  slow  variable,  the  stored  CRH,  authors  were  able  to 
observe  ultradian  oscillations  in  cortisol  level  and  stress- 
induced  transitions  into  a  low-cortisol  diseased  state.  The 
manuscript  is  clearly  written  and  I  recommend  a  publica¬ 
tion  in  Biology  Direct. 

Authors'  response:  We  thank  Dr.  Ha  Youn  Lee  for  her 
positive  recommendation. 

Recommendations:  The  validity  of  this  model  can  be 
more  clearly  demonstrated  by  comparing  the  model  solu¬ 
tion  with  the  previously  published  data  of  cortisol  dynam¬ 
ics  in  normal  and  post-traumatic  stress  disorder  (PTSD) 
in  references  ([53]  and  [16]). 

Authors'  response:  We  thank  Dr.  Lee  for  the  sugges¬ 
tion.  It  would  indeed  be  interesting  to  fit  our  model  to 
the  data  in  [S3]  and  compare  the  estimated  parameter 
values  to  those  of  Siriam  et  al.  [16],  as  their  estima¬ 
tions  were  also  based  on  the  data  from  [S3].  As  noted 
in  the  manuscript ,  Siriam  et  al.  have  estimated  Ki,  the 
parameter  analogous  to  I<a  in  our  model  that  represents 
the  strength  of  negative  feedback  imparted  by  cortisol 
in  the  pituitary.  Their  estimations  of  I<i  were  consistent 


with  an  enhanced  negative  feedback  action  of  cortisol  (i.e. 
decreased  Ki)  in  PTSD  patients,  as  hypothesized  by  Yehuda 
etal.  [S3]. 

Instead  of  fitting  our  model  to  the  limited  number  of  data 
( time  series  data  of  three  individuals,  one  from  each  of 
control,  depressed,  and  PTSD  group)  used  in  [16],  we  can 
predict  how  our  model  solution  will  compare  to  the  results 
of  [16]  by  analyzing  the  effect  of  varying  the  analogous 
parameter  (Ka)  on  the  nucllcine  structure  of  our  model. 
When  Ka  is  decreased  (enhanced  negative  feedback  in  the 
pituitary),  the  q-nullcline  in  our  model  shifts  toward  the 
normal,  high  cortisol  state  (on  the  upper  branch  of  the  c- 
nullcline),  away  from  the  diseased  state.  Thus,  contrary 
to  the  hypothesis  of  Yehuda  et  al.  [S3],  enhanced  negative 
feedback  action  of  cortisol  does  not  characterize  low  cor¬ 
tisol  levels  observed  in  PTSD  patients  within  our  model. 
We  plan  to  provide  further  discussion  on  this  matter  in  our 
future  work. 

Minor  issues:  The  observation  that  transition  to  dis¬ 
eased  state  depends  on  stress  timing  is  interesting  but 
it  can  be  addressed  whether  or  not  this  is  a  biologically 
relevant  phenomenon. 

Authors'  response:  In  a  previous  study  [21  ],  it  was  shown 
that  changes  in  corticosterone  levels  induced  by  acute  audi¬ 
tory  stressor  in  rats  were  dependent  on  the  timing  of  stress 
onset  relative  to  the  phase  of  underlying  corticosterone 
oscillations.  Based  on  this  observation,  we  believe  that  the 
timing  of  stress  onset  could  be  relevant  in  the  transitions 
to  diseased  states.  Please  refer  to  the  Additional  file  1 
for  details  on  how  the  experimental  observation  can  be 
explained  within  our  model. 

To  address  the  timing  issue  in  a  more  detailed  manner, 
we  need  a  better  description  of  the  synaptic  input  func¬ 
tion  of  the  PVN,  I(t),  that  models  how  and  when  stress 
response  is  initiated  and  terminated.  We  have  shown  in 
the  manuscript  that  timing  may  be  crucial  in  inducing 
transitions,  but  only  when  the  stress  duration  is  near  its 
transition  thresholds.  A  more  realistic  form  of  I(t)  will 
thus  allow  us  to  understand  under  what  circumstances 
stress  duration  may  be  near  such  transition  thresholds.  We 
are  are  currently  investigating  the  endocannabinoid  system 
that  is  known  to  regulate  the  initiation  and  termination  of 
stress  response,  which  is  subject  to  fast  nongenomic  actions 
of  cortisol. 

Additional  file 


Additional  file  1 :  The  file  provides  appendices  describing  the 
mathematical  details  of  our  analysis.  It  includes  discussions  on 
nondimensionalization  and  parameter  estimations  used  in  our  model. 
Analyses  of  the  nullcline  structure  and  the  parameter  space  are  also 
provided  and  the  effect  of  stress  onset  timing  on  initial  stress  response  is 
described  in  detail.  Finally,  a  possible  form  of  cortisol  dependent  /ex t  is 
illustrated  with  its  implications  on  model  behavior.  (PDF  2099  kb) 
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Abbreviations 

ACTH:  adrenocorticotropic  hormone;  BEP:  brief  eclectic  psychotherapy;  CFS: 
chronic  fatigue  syndrome;  CORT:  glucocorticoids;  CRH:  corticotropin-releasing 
hormone;  CRHR-1 :  Corticotropin-releasing  hormone  receptor  1 ;  DDE:  delay 
differential  equation;  ET:  exposure  therapy;  GH:  growth  hormone;  GnRH: 
gonadotropin-releasing  hormone;  GR:  glucocorticoid  receptors;  MDD:  major 
depressive  disorder;  hnRNA:  heterogeneous  nuclear  RNA;  HPA: 
hypothalamic-pituitary-adrenal;  HPG:  hypothalamic-pituitary-gonadal;  mRNA: 
messenger  ribonucleic  acid;  ODE:  ordinary  differential  equation;  PA: 
pituitary-adrenal;  POA:  preoptic  area;  PTSD:  post-traumatic  stress  disorder; 

PVN:  paraventricular  nucleus;  TBI:  traumatic  brain  injury. 

Competing  interests 

The  authors  declare  that  they  have  no  competing  interests. 

Authors'  contributions 

LK  performed  the  mathematical  analysis  and  computational  simulations.  MRD 
and  TC  formulated  and  designed  the  study.  All  authors  contributed 
significantly  to  the  analysis  of  the  mathematical  simulations  and  the  writing  of 
the  manuscript.  All  authors  read  and  approved  the  final  manuscript. 

Acknowledgments 

This  work  was  supported  by  the  Army  Research  Office  via  grant 

W91 1 NF-14-1  -0472  and  the  NSF  through  grant  BCS-1 3481 23.  The  authors  also 

thank  professors  T.  Minor  and  M.  Wechselberger  for  insightful  discussions. 

Author  details 

1  Department  of  Biomathematics,  Univ  of  California,  Los  Angeles,  51 09  Life 
Sciences  621  Charles  E.  Young  Dr.  South,  Los  Angeles,  USA.  department  of 
Mathematics,  CalState-Northridge,  18111  Nordhoff  St.,  Los  Angeles,  USA. 

3 Department  of  Biomathematics  and  Department  of  Mathematics,  University 
of  California,  Los  Angeles,  5209  Life  Sciences  621  Charles  E.  Young  Dr.  South, 
Los  Angeles,  USA. 

Received:  24  December  201 5  Accepted:  1 7  March  201 6 
Published  online:  25  March  2016 

References 

1 .  Denver  R.  Structural  and  functional  evolution  of  vertebrate 
neuroendocrine  stress  systems.  Ann  N  Y  Acad  Sci.  2009;1 163(1  ):1  -16. 

2.  Gold  P,  Chrousos  G.  Organization  of  the  stress  system  and  its 
dysregulation  in  melancholic  and  atypical  depression:  high  vs  low 
CRH/NE  states.  Mol  Psychiatry.  2002;7(3):254-75. 

3.  Juruena  M,  CleareA,  Pariante  C.  The  hypothalamic  pituitary  adrenal  axis, 
glucocorticoid  receptor  function  and  relevance  to  depression.  Rev  Bras 
Psiquiatr.  2004;26(3):1 89-201 . 

4.  RohlederN,  Joksimovic  L,  WolfJ,  Kirschbaum  C.  Hypocortisolism  and 
increased  glucocorticoid  sensitivity  of  pro-inflammatory  cytokine 
production  in  bosnian  war  refugees  with  posttraumatic  stress  disorder. 
Biol  Psychiatry.  2004;55(7):745-51. 

5.  Giorgio  AD,  Hudson  M,  JerjesW,  Cleare  A.  24-hour  pituitary  and  adrenal 
hormone  profiles  in  chronic  fatigue  syndrome.  Psychosom  Med. 
2005;67(3):433-40. 

6.  Jerjesnd  W,  Peters  T,  Taylor  N,  Wood  P,  WesselyS,  Cleare  A.  Diurnal 
excretion  of  urinary  cortisol,  cortisone,  and  cortisol  metabolites  in  chronic 
fatigue  syndrome.  J  Psychosom  Res.  2006;60(2):  145-53. 

7.  Crofford  L,  Young  E,  Cary  NEK,  Korszun  A,  Brucksch  C,  McClure  L, 

Brown  M,  Demitrack  M.  Basal  circadian  and  pulsatile  ACTH  and  cortisol 
secretion  in  patients  with  fibromyalgia  and/or  chronic  fatigue  syndrome. 
Brain  Behav  Immun.  2004;18(4):3 14-25. 

8.  Yehuda  R,  TeicherM,  Levengood  R,  Trestman  R,  Siever  L.  Circadian 
regulation  of  basal  cortisol  levels  in  posttraumatic  stress  disorder.  Ann  N  Y 
Acad  Sci.  1994;746(1):378-80. 

9.  Vinther  F,  Andersen  M,  Ottesen  JT.  The  minimal  model  of  the 
hypothalamic-pituitary-adrenal  axis.  J  Math  Biol.  201 1  ;63:663— 90. 

10.  JelicS,  CupicZ,  Kolar-Anic  L.  Mathematical  modeling  of  the 
hypothalmic-pituitary-adrenal  system  activity.  Math  Biosci.  2005;197: 
173-87. 

1 1 .  Kyrylov  V,  Severyanova  L,  Vieira  A.  Modeling  robust  oscillatory  behavior 
of  the  hypothalamic-pituitary-adrenal  axis.  IEEE  Trans  Biomed  Eng. 
2005;52(12):1 977-83. 


1 2.  Savic  D,  Knezevic  G,  Opacic  G.  A  mathematical  model  of  stress  reaction: 
Individual  differences  in  threshold  and  duration.  Psychobiology. 
2000;28(4):581-92. 

13.  Walker  JJ,  Terry  JR,  Lightman  SL.  Origin  of  ultradian  pulsatility  in  the 
hypothalamic-pituitary-adrenal  axis.  Proc  R  Soc  Lond  B  Biol  Sci. 

201 0;277(1 688):  1 627-33. 

14.  Rankin  J,  Walker  J,  WindleR,  Lightman  S,  Terry  J.  Characterizing  dynamic 
interactions  between  ultradian  glucocorticoid  rhythmicity  and  acute 
stress  using  the  phase  response  curve.  PloS  ONE.  201 2;7(2):30978. 

1 5.  Bairagi  N,  Chatterjee  S,  Chattopadhyay  J.  Variability  in  the  secretion  of 
corticotropin-releasing  hormone,  adrenocorticotropic  hormone  and 
cortisol  and  understandability  of  the  hypothalamic-pituitary-adrenal  axis 
dynamics  —  a  mathematical  study  based  on  clinical  evidence.  Math  Med 
Biol.  2008;25:37-63. 

16.  Sriram  K,  Rodriguez-Fernandez  M,  Doyle  III  FJ.  Modeling  cortisol 
dynamics  in  the  neuro-endocrine  axis  distinguishes  normal,  depression, 
and  post-traumatic  stress  disorder  (PTSD)  in  humans.  PLoS  Comput  Biol. 
201 2;8:1 002379. 

17.  Gupta  S,  Aslakson  E,  Gurbaxani  BM,  Vernon  SD.  Inclusion  of  the 
glucocorticoid  receptor  in  a  hypothalamic  pituitary  adrenal  axis  model 
reveals  bistability.  Theor  Biol  Med  Model.  2007;4:8. 

18.  WindleR,  Woods,  Lightman  S,  Ingram  C.  The  pulsatile  characteristics  of 
hypothalamo-pituitary-adrenal  activity  in  female  Lewis  and  Fischer  344 
rats  and  its  relationship  to  differential  stress  responses.  Endocrinology. 

1 998;1 39(1 0):4044-52. 

1 9.  Chrousos  G.  Editorial:  ultradian,  circadian,  and  stress-related 
hypothalamic-pituitary-adrenal  axis  activity  —  a  dynamic 
digital-to-analog  modulation.  Endocrinology.  1998;139(2):437-40. 

20.  Conway-Campbell  B,  Sarabdjitsingh  R,  McKenna  M,  PooleyJ,  Kershaw  Y, 
MeijerO,  Kloet  ED,  Lightman  S.  Glucocorticoid  ultradian  rhythmicity 
directs  cyclical  gene  pulsing  of  the  clock  gene  period  1  in  rat 
hippocampus.  J  Neuroendocrinol.  201 0;22(1 0):  1 093—1 00. 

21.  WindleR,  Woods,  Shanks  N,  Lightman  S,  Ingram  C.  Ultradian  rhythm  of 
basal  corticosterone  release  in  the  female  rat:  Dynamic  interaction  with 
the  response  to  acute  stress.  Endocrinology.  1 998;1 39(2):443-50. 

22.  Watts  A.  Glucocorticoid  regulation  of  peptide  genes  in  neuroendocrine 
CRH  neurons:  a  complexity  beyond  negative  feedback.  Front 
Neuroendocrinol.  2005;26(3):1 09-30. 

23.  Ono  N,  Castro  JD,  McCann  S.  Ultrashort-loop  positive  feedback  of 
corticotropin  (ACTH)-releasing  factor  to  enhance  ACTH  release  in  stress. 
Proc  Natl  Acad  Sci.  1 985;82(1 0):3528-31 . 

24.  FitzHugh  R.  Mathematical  models  of  threshold  phenomena  in  the  nerve 
membrane.  Bull  Math  Biophys.  1955;17(4):257-78. 

25.  Silva  FLD,  BlanesW,  KalitzinS,  Parra  J,  Suffczynski  P,  Velis  D.  Epilepsies  as 
dynamical  diseases  of  brain  systems:  basic  models  of  the  transition 
between  normal  and  epileptic  activity.  Epilepsia.  2003;44(s1 2):72-83. 

26.  Ben-Zvi  A,  Vernon  SD,  Broderick  G.  Model-based  therapeutic  correction 
of  hypothalamic-pituitary-adrenal  axis  dysfunction.  PLoS  Comput  Biol. 
2009;5(1):1 000273. 

27.  Tsai  SY,  Carlstedt-Duke  J,  Weigel  NL,  Dahlman  K,  Gustafsson  J,  Tsai  MJ, 
O'Malley  BW.  Molecular  interactions  of  steroid  hormone  receptor  with  its 
enhancer  element:  evidence  for  receptor  dimer  formation.  Cell. 

1 988;55(2):36 1  —9. 

28.  Andersen  M,  Vinther  F,  Ottesen  J.  Mathematical  modeling  of  the 
hypothalamic-pituitary-adrenal  gland  (HPA)  axis,  including  hippocampal 
mechanisms.  Math  Biosci.  201 3;246(1  ):1 22-38. 

29.  Papaikonomou  E.  Rat  adrenocortical  dynamics.  J  Physiol.  1 977;265(1 ): 
119-31. 

30.  EnglerD,  PhamT,  LiuJ,  Fullerton  M,  Clarke  I,  Funder  J.  Studies  of  the 
regulation  of  the  hypothalamic-pituitary-adrenal  axis  in  sheep  with 
hypothalamic-pituitary  disconnection.  II,  evidence  for  in  vivo  ultradian 
hypersecretion  of  proopiomelanocortin  peptides  by  the  isolated 
anterior  and  intermediate  pituitary.  Endocrinology.  1 990;1 27(4): 

1956-66. 

31.  WeiserM,  OsterlundC,  Spencer  R.  Inhibitory  effects  of  corticosterone  in 
the  hypothalamic  paraventricular  nucleus  (pvn)  on  stress-induced 
adrenocorticotrophic  hormone  secretion  and  gene  expression  in  the  pvn 
and  anterior  pituitary.  J  Neuroendocrinol.  201 1;23(1 2):1 231-40. 

32.  Ma  X,  Aguilera  G.  Differential  regulation  of  corticotropin-releasing 
hormone  and  vasopressin  transcription  by  glucocorticoids. 
Endocrinology.  1 999;1 40(1 2):5642-50. 


Ki  m  et  al.  Biology  Direct  (20 1 6)  1 1 : 1 3 


Page  1 8  of  1 8 


33.  Watts  A,  Sanchez-Watts  G.  Region-specific  regulation  of  neuropeptide 
mRNAs  in  rat  limbic  forebrain  neurones  by  aldosterone  and 
corticosterone.  J  Physiol.  1995;484(3):721-36. 

34.  Tasker  J,  Di  S,  Malcher-Lopes  R.  Rapid  glucocorticoid  signaling  via 
membrane-associated  receptors.  Endocrinology.  2006;1 47(1 2): 

5549-56. 

35.  Kasai  M,  Yamashita  H.  Inhibition  by  cortisol  of  neurons  in  the 
paraventricular  nucleus  of  the  hypothalamus  in  adrenalectomized  rats;  an 
in  vitro  study.  Neurosci  Lett.  1 988;91  (1  ):59— 64. 

36.  Kasai  M,  Yamashita  H.  Cortisol  suppresses  noradrenaline-induced 
excitatory  responses  of  neurons  in  the  paraventricular  nucleus;  an  in  vitro 
study.  Neurosci  Lett.  1 988;91  (1  ):65— 70. 

37.  Jones  M,  HillhouseE,  Burden  J.  Dynamics  and  mechanics  of 
corticosteroid  feedback  at  the  hypothalamus  and  anterior  pituitary  gland. 
J  Endocrinol.  1 977;73(3):405-1 7. 

38.  Ginsberg  A,  Campeau  S,  Day  H,  Spencer  R.  Acute  glucocorticoid 
pretreatment  suppresses  stress-induced  hypothalamic-pituitary-adrenal 
axis  hormone  secretion  and  expression  of  corticotropin-releasing 
hormone  hnRNA  but  does  not  affect  c-fos  mRNA  or  fos  protein 
expression  in  the  paraventricular  nucleus  of  the  hypothalamus.  J 
Neuroendocrinol.  2003;1 5(1  1  ):1 075-83. 

39.  Chen  Y,  Hua  S,  Wang  C,  Wu  L,  Gu  Q,  Xing  B.  An  electrophysiological 
study  on  the  membrane  receptor-mediated  action  of  glucocorticoids  in 
mammalian  neurons.  Neuroendocrinology.  1 991  ;53(Suppl.  1  ):25— 30. 

40.  ImakiT,  Xiao-Quan  W,  ShibasakiT,  Yamada  K,  Harada  S,  Chikada  N, 
Naruse  M,  Demura  H.  Stress-induced  activation  of  neuronal  activity  and 
corticotropin-releasing  factor  gene  expression  in  the  paraventricular 
nucleus  is  modulated  by  glucocorticoids  in  rats.  J  Clin  Investig.  1 995;96(1 ): 
231. 

41.  Papadimitriou  A,  Priftis  K.  Regulation  of  the 
hypothalamic-pituitary-adrenal  axis.  Neuroimmunomodulation. 

2009;1 6(5):265. 

42.  MakinoS,  Hashimoto  K,  Gold  P.  Multiple  feedback  mechanisms 
activating  corticotropin-releasing  hormone  system  in  the  brain  during 
stress.  Pharmacol  Biochem  Behav.  2002;73(1):1 47-58. 

43.  Herman  JP,  Figueiredo  H,  Mueller  NK,  Ulrich-LaiY,  Ostrander  M,  Choi  D, 
Cullinan  W.  Central  mechanisms  of  stress  integration:  hierarchical  circuitry 
controlling  hypothalamo-pituitary-adrenocortical  responsiveness.  Front 
Neuroendocrinol.  2003;24(3):  151  —80. 

44.  McEwen  BS,  Stellar  E.  Stress  and  the  individual:  mechanisms  leading  to 
disease.  Arch  Intern  Med.  1 993;1 53:2093-1 01 . 

45.  McEwen  BS.  Stress,  adaptation,  and  disease:  Allostasis  and  allostatic  load. 
Ann  N  Y  Acad  Sci.  1 998;840(1  ):33-44. 

46.  Dince  SM,  Rome  RD,  McEwen  BS,  Tang  AC.  Enhancing  offspring 
hypothalamic-pituitary-adrenal  (hpa)  regulation  via  systematic  novelty 
exposure:  the  influence  of  maternal  HPA  function.  Front  Behav  Neurosci. 
201 4;8:204.doi:10.3389/fnbeh.201 4.00204. 

47.  Hauger  RL,  Thrivikraman  KV,  Plotsky  PM.  Age-related  alterations  of 
hypothalamic-pituitary-adrenal  axis  function  in  male  Fischer  344  rats. 
Endocrinology.  1 994;1 34(3):  1 528-36. 

48.  Averill  P,  Beck  J.  Posttraumatic  stress  disorder  in  older  adults:  a 
conceptual  review.  J  Anxiety  Disord.  2000;1 4(2):  133-56. 

49.  RegierD,  BoydJ,  Burke  J,  Rae  D,  Myers  J,  Kramer  M,  Robins  L,  George  L, 
Kamo  M,  Locke  B.  One-month  prevalence  of  mental  disorders  in  the 
united  states:  based  on  five  epidemiologic  catchment  area  sites.  Arch 
Gen  Psychiatr.  1 988;45(1 1  ):977-86. 

50.  Simerly  RB,  Swanson  LW,  Chang  C,  Muramatsu  M.  Distribution  of 
androgen  and  estrogen  receptor  mrna-containing  cells  in  the  rat  brain: 

An  in  situ  hybridization  study.  J  Comp  Neurol.  1 990;294(1  ):76— 95. 

51.  Greco  B,  Allegretto  E,  Tetel  M,  Blaustein  J.  Coexpression  of  ER ft  with  ERa 
and  progestin  receptor  proteins  in  the  female  rat  forebrain:  effects  of 
estradiol  treatment.  Endocrinology.  2001  ;1 42(1 2):51 72-81. 

52.  Feldman  S,  Conforti  N,  Saphier  D.  The  preoptic  area  and  bed  nucleus  of 
the  stria  terminalis  are  involved  in  the  effects  of  the  amygdala  on 
adrenocortical  secretion.  Neuroscience.  1 990;37(3):775— 9. 

53.  Yehuda  R,  TeicherM,  Levengood  R,  Trestman  R,  Levengood  R,  Siever  L. 
Cortisol  regulation  in  posttraumatic  stress  disorder  and  major  depression: 
a  chronobiological  analysis.  Biol  Psychiatry.  1996;40(2):79-88. 

54.  Yehuda  R,  LeDoux  J.  Response  variation  following  trauma:  a  translational 
neuroscience  approach  to  understanding  PTSD.  Neuron.  2007;56: 

1 9-32. 


55.  OlffM,  de  Vries  G,  GuzelcanY,  Assies  J,  Gersons  B.  Changes  in  cortisol 
and  DHEA  plasma  levels  after  psychotherapy  for  PTSD. 
Psychoneuroendocrinology.  2007;32(6):61 9-26. 

56.  Foa  E,  Keane  T,  Friedman  M,  Cohen  J.  Effective  treatments  for  PTSD: 
practice  guidelines  from  the  international  society  for  traumatic  stress 
studies.  New  York:  Guilford  Press;  2008. 

57.  Rauch  S,  Eftekhari  A,  Ruzek  J.  Review  of  exposure  therapy:  a  gold 
standard  for  PTSD  treatment.  J  Rehabil  Res  Dev.  201 2;49(5):679-88. 

58.  Trouche  S,  Sasaki  J,  Tu  T,  Reijmers  L.  Fear  extinction  causes 
target-specific  remodeling  of  perisomatic  inhibitory  synapses.  Neuron. 
201 3;80(4):1 054-65. 

59.  Lemieux  A,  Coe  C.  Abuse-related  posttraumatic  stress  disorder:  evidence 
for  chronic  neuroendocrine  activation  in  women.  Psychosom  Med. 

1 995;57(2):1 05-1 5. 

60.  Baum  A.  Implications  of  psychological  research  on  stress  and 
technological  accidents.  Am  Psychol.  1 993;48(6):665. 

61 .  Rosenmund  C,  Stevens  C.  Definition  of  the  readily  releasable  pool  of 
vesicles  at  hippocampal  synapses.  Neuron.  1 996;1 6(6):1 1 97-207. 


Submit  your  next  manuscript  to  BioMed  Central 
and  we  will  help  you  at  every  step: 

•  We  accept  pre-submission  inquiries 

•  Our  selector  tool  helps  you  to  find  the  most  relevant  journal 

•  We  provide  round  the  clock  customer  support 

•  Convenient  online  submission 

•  Thorough  peer  review 

•  Inclusion  in  PubMed  and  all  major  indexing  services 

•  Maximum  visibility  for  your  research 

Submit  your  manuscript  at  ^  s 

www.biomedcentral.com/submit  BioMed  Central 


